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20.  Abstract 


Gas-phase  chemical  reactions  play  an  important  role  in  determining  the 
power  and  sensitivity  of  energetic  materials  which  are  used  as  fuels 
in  military  propulsion  systems.  Theoretical  methods  can  provide  detailed 
dynamical  information  about  the  important  chemical  reactions,  thereby 
aiding  in  the  design  of  new  propellants.  High  order  electronic  structure 
calculations  with  empirical  bond  additivity  corrections  have  been  applied 
to  prototype  reactions  and  validated  by  comparisons  of  different  levels 
of  theory  and  by  comparison  with  experiment.  These  methods  provide  the 
most  cost  effective  means  of  obtaining  thermal  rate  data  and  show  promise 
of  being  applicable  to  reactions  involving  large  polyatomic  molecules 
important  in  nitrocubane  decomposition. 
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EXECUTIVE  SUMMARY 


This  document  reports  the  results  of  a  six-month  Phase  I  small  business  innovative 
research  program  through  the  strategic  defense  initiative.  The  main  objective  of  this  program  was 
to  establish  the  feasibility  of  using  theoretical  methods  for  obtaining  kinetic  data  important  in 
nitrocubane  decomposition.  The  potential  utility  of  theoretical  methods  for  the  aid  in  the  design  of 
new  propellants  is  unquestionable;  what  remains  in  question  is  whether  the  methods  are  of 
sufficient  accuracy  and  capable  of  treating  large  enough  molecules  to  be  of  use.  The  present  report 
outlines  the  first  successful  steps  taken  in  answering  these  questions. 

The  research  preformed  on  this  program  entailed  a  collaborative  effort  of  electronic 
structure  and  dynamics  calculations.  The  electronic  structure  calculations  provide  information 
about  the  potential  energy  surface  which  is  then  used  in  dynamics  calculations  of  the  rate  constants. 
The  sensitivity  of  the  computed  rate  constants  to  the  level  of  theory  was  studied  for  a  relatively 
small  reaction  (involving  five  atoms)  for  which  high  levels  of  theory  could  be  afforded. 
Comparison  was  also  made  between  the  computed  and  experimental  rates  as  further  validation  of 
the  methods.  These  studies  indicate  that  accurate  rate  constants  (to  within  ab out  50%)  can  be 
obtained  for  a  modest  computational  effort  -  several  hours  of  Cray  Is  computer  time  were  required 
to  obtain  rates  for  temperatures  from  200  to  2400  K. 

In  order  to  treat  reactions  involving  very  large  molecules  (such  as  the  initial  unimolecular 
decomposition  step  in  nitrocubane  combustion)  approximate  methods  for  estimating  kinetic 
parameters  are  needed.  As  a  first  step  in  this  direction  we  have  carried  out  a  systematic  study  of 
hydrogen  attack  on  a  series  of  nitro  containing  molecules  to  examine  the  effects  of  substituents  on 
the  computed  reaction  rates.  In  addition,  a  new  thermochemical  kinetic  analysis  has  been  presented 
which  includes  important  dynamical  effects  that  are  cmcial  for  reliable  predictions  of  thermal  rates. 
This  new  analysis  will  be  the  basis  for  methods  of  estimating  kinetic  parameters  for  large 
polyatomic  systems  from  bond  and  group  additivity  relationships,  but  using  information  about 
potential  enemy  surfaces  for  smaller  analogous  reactions. 
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Appendices 

A.  Theoretical  calculations  of  the  thermal  rate  constants  for  the  gas-phase 
chemical  reactions  11  +  NH-j  ^  H2  +  NH2  and  D  +  ND3  — »  D2  +  ND2 

B .  Thermochemical  kinetic  analysis  of  tunneling  and  the  incorporation  of 
tunneling  contributions  in  thermochemical  kinetics 


1.  Introduction 


Propulsion  requirements  for  the  strategic  defense  initiative  require  the  development  of  new 
fuels  that  are  lightweight  yet  highly  energetic.  Fuels  based  upon  combining  beryllium  with  an 
oxidizer  show  promise  of  meeting  this  criteria  of  a  high  specific  impulse.  The  ideal  oxidizing  agent 
is  one  which  is  highly  energetic  yet  stable  and,  therefore,  detonates  in  a  well  characterized, 
predictable  manner.  An  example  of  a  highly  energetic  yet  stable  molecule  is  the  cage  compound 
cubane1 


with  much  of  its  energy  arising  from  the  ring  strain  imposed  by  the  cubic  geometry  of  the  carbon 
skeleton.  There  is  currently  a  large  effort  within  SDI  to  develop  high-density  oxidizers  by 
substituting  cubanes  with  nitro  (NO-,)  or  nitro-containing  groups.  The  explosive  nature  of 

nitrocubanes  derives  from  the  ring  strain  and  the  energy-rich  oxidizing  nitro  substituents. 
Nitrocubanes  containing  many  nitro  groups  are  more  energetic  and  more  effective  oxidizers,  but 
also  tend  to  to  be  more  unstable.  Some  stability  can  be  regained  by  placing  the  nitro  groups  onto 
substituents  that  move  them  further  away  from  the  cube.  Examples  of  recent  compounds  along 
this  line  are  the  ammonium  nitronate  salts,  e.g., 

NH3C(N02)3 

The  synthesis  and  characterization  of  nitrocubanes  such  as  these  is  currently  of  great  interest  and  is 
being  pursued  in  several  laboratories. 

One  of  the  factors  controlling  the  stability  of  nitrocubanes  is  the  rate  of  the  unimolecular 
decomposition  of  the  parent  molecule  in  the  condensed  phase.  This  is  controlled  largely  by  the 
bond  strength  of  the  weakest  bond  in  the  molecule.  To  fully  realize  the  power  stored  in  these 
molecules  they  must  be  converted  into  more  thermodynamically  stable  molecules  such  as  H-,0. 


CO-,,  and  N->.  This  conversion  is  controlled  by  chemistry  occurring  in  the  gas  phase.  Gas-phase 
chemical  reactions  play  an  important  role  in  the  sensitivity  as  well  as  the  energetics  of  the 
propellants.  The  exothermicity  of  the  gas-phase  reactions  determine  their  energy  output  and  the 
kinetics  of  the  gas-phase  reactions  control  the  bum  rate,  and  together  they  determine  the  power  of 
the  propellant.2  Also,  radicals  which  are  products  of  the  gas-phase  reactions  can  attack  the  parent 
molecules  and  affect  the  sensitivity  of  the  nitrocubane.  The  effect  of  substituent  groups  on  the 
sensitivity  of  nitrocubanes  is  twofold.  Firstly,  the  substituent  can  affect  the  sensitivity  by  changing 
the  strength  of  the  weakest  bond  in  the  parent  molecule,  and  secondly,  the  substituent  can  control 
the  sensitivity  by  altering  the  mix  of  H/C/N/O  which  affects  the  gas-phase  kinetics.  An 
understanding  of  the  kinetically  important  reaction  pathways  which  are  operative  in  the  complex 
decomposition  mechanism  is  crucial  to  developing  energetic  yet  stable  oxidizers. 

Theoretical  methods  can  make  a  contribution  to  the  understanding  of  the  decomposition  of 
propellants.  To  date,  great  strides  have  been  made  in  computer  modeling  of  the  ignition  and  flames 
of  propellants  such  as  RDX.3  Studies  of  ignition  solve  the  nonlinear  coupled  differential  equations 
for  temperature  and  species  mass  fraction  as  a  function  of  time.  Comparison  with  experimental 
observations  of  gross  features  of  the  ignition  (such  as  ignition  time  or  observed  formation  rates  of 
specific  molecules)  tests  the  underlying  reaction  mechanism  and  kinetic  data  used  as  input  into  the 
computer  model.  The  level  of  detailed  species  concentration  as  a  function  of  time  that  is  afforded 
by  these  simulations  is  currently  not  available  in  experimental  studies.  This  detail  enables 
identification  of  important  reaction  pathways  -  the  important  elementary  reaction  steps  -  and 
sensitivity  analysis  identifies  crucial  kinetic  data  that  must  be  determined  accurately. 

Great  strides  have  also  been  made  in  the  theoretical  methods  for  describing  the  kinetics  of 
elementary  reactions.  We  are  currently  on  the  threshold  of  attaining  the  necessary  capabilities  to 
enable  routine  calculations  of  gas-phase  chemical  reaction  rates.  With  this  capability  in  hand,  it 
will  be  easy  to  test  the  influence  of  substituents  on  the  bond  strengths  and,  thereby,  unimolecular 
deeomposition  rates,  and  also  to  test  their  influence  on  the  subsequent  gas-phase  kinetics. 
Theoretical  calculations  provide  an  important  complementary  tool  to  experimental  gas  kinetic 


studios.  Elementary  reaction  steps  can  he  studied  theoretically  for  reactant  species  that  are  very 
short  lived  and  hard  to  prepare  experimentally.  Furthermore,  the  theoretical  calculations  can  be 
toed  to  extend  experimental  data  to  physical  conditions  which  are  not  attainable  in  the  laboratory. 

In  collaboration  with  computer  modeling  efforts,  these  theoretical  studies  can  provide  physical 
insight  into  the  reaction  pathways  and  detailed  dynamics  of  the  chemical  reactions.  For  systems 
w  hich  are  not  too  large  the  theoretical  methods  are  capable  of  quantitative  predictions  of  the  gas- 
phase  reaction  rates  and  for  very  large  systems  these  methods  will  give  qualitative  rate  data  which 
w  ill  allow  the  assessment  of  the  most  important  reaction  pathways.  These  studies  are  crucial  to  the 
detailed  understanding  of  substituent  effects  upon  the  energetics  and  sensitivity  of  the  oxidizers. 

One  objective  of  the  Phase  I  research  has  been  to  extend  and  validate  the  theoretical 
methods  used  for  computing  rate  constants  for  elementary  gas-phase  reactions  and  to  demonstrate 
the  feasibility  of  a  theoretical  approach  for  describing  the  dynamics  of  these  gas-phase  reactions. 
The  theoretical  calculation  of  the  chemical  reaction  rate  is  a  two  step  process:  first,  the  interaction 
energies  between  the  atoms  are  obtained  from  an  electronic  structure  (quantum  chemistry) 
calculation,  followed  by  a  dynamical  calculation  of  the  rate.  Although  the  electronic  structure 
calculations  do  provide  vital  insight  into  the  energetically  accessible  reaction  pathways,  they 
provide  only  structural  information  which  is  often  not  enough;  dynamical  constraints  can  make 
energetically  accessible  reaction  pathways  so  slow  that  they  become  unimportant  in  the  kinetic 
mechanism.  A  collaborative  effort  which  successfully  interfaces  the  quantum  chemistry  and 
dynamics  calculations  is  needed  to  accurately  determine  the  dynamics  of  elementary  gas-phase 
reactions. 

The  d\  namies  cal.  alations  are  based  upon  variational  transition  state  theory  with  quantum 
mechanical  tunneling  effects  included  by  semidassical  adiabatic  ground-state  transmission 
oxM'ficients.4  I'he  use  of  variational  transition  slate  theory  in  these  calculations  is  mandated  for 
several  reasons:  <  1 )  it  includes  the  factors  most  important  in  controlling  the  rate  of  chemical 
reactions  and  is  currently  the  most  cost  effective  method  for  obtaining  reliable  predictions  of  rates 
for  a  variety  of  gas-phase  reactions;  (2)  it  is  capable  of  utilizing  limited  information  about  the 
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potential  along  the  minimum  energy  path  without  requiring  a  global  potential  energy  surface;  (3) 
because  of  the  nature  of  the  semiempirical  potential  information  utilized  in  this  work,  the  dynamical 
bottleneck  for  the  reaction  does  not  occur  at  the  saddle  point  for  the  reaction  and  a  variational 
procedure  to  locate  the  dynamical  bottleneck  is  needed;  and  (4)  it  provides  a  consistent  method  for 
incorporating  quantum  mechanical  tunneling  effects  which  are  crucial  for  accurate  predictions  of 
the  rates,  especially  for  reactions  which  actively  involve  hydrogen  atoms  at  temperatures  below 
600  K. 

Rigorous  ab  initio  electronic  structure  calculations  are  currently  incapable  of  predicting  the 
energetics  of  a  polyatomic  reaction  accurately  enough  to  allow  even  a  qualitative  estimate  of  the  rate 
constant.  The  most  successful  electronic  structure  methods  incorporate  some  type  of  correction, 
justified  either  on  a  theoretical  or  empirical  basis.  One  such  approach  is  the  BAC-MP4  method5 
developed  by  Dr.  Carl  Melius  at  Sandia  National  Laboratories,  Livermore,  California  (SNLL). 

The  present  research  was  performed  in  collaboration  with  Dr.  Melius  and  Dr.  M.  L.  Koszykowski 
at  SNLL.  The  BAC-MP4  method  has  been  successfully  applied  to  the  calculation  of  the 
thermochemistry  of  bound  species  and  is  currently  the  foremost  method  for  prediction  of 
thermochemical  data  for  transient  species.  In  the  present  application  the  BAC-MP4  method  is 
extended  to  calculations  of  potential  information  in  the  interaction  region  of  the  potential  rather  than 
at  equilibrium  geometries.  A  major  goal  of  this  work  is  the  validation  of  this  method  for  obtaining 
the  potential  information  needed  for  the  dynamics  calculations. 

Another  objective  of  this  research  was  to  treat  chemical  reactions  of  current  interest  and 
importance  in  propellant  chemistry.  In  this  direction  it  was  proposed  to  validate  the  theoretical 
methods  in  a  study  of  the  simplest  prototype  of  a  radical  attack  on  a  nitro  group 
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I 


R  +  HONO 


(R2) 


RNO  +  OH 


(R3) 


where  R  is  replace  by  a  hydrogen  atom.  Towards  this  end,  we  examined  the  prototypical  reactions 
H  +  NH, ->  H2  +  NH2  (R4) 

H,  +  NH2 H  +  NH3  (R5) 

and 

D  +  ND-,  — >  Dt  +  ND-,  (R6) 

These  reactions  provide  a  wealth  of  experimental  data  for  a  critical  test  of  the  theoretical  methods 
and  the  reactions  are  small  enough  (have  a  small  enough  number  of  electrons)  to  afford  relatively 
large  electronic  structure  calculations  for  the  validation.  Subsequently  reactions  (R1)-(R3)  were 
studied  for  R  =  H,  CH3,  OH,  NHi,  and  CH3NH  but  no  validation  of  the  theoretical  methods  was 

attempted  for  these  reactions. 

The  Phase  I  objectives  have  been  surpassed  in  the  six  month  research  effort.  The 
feasibility  of  the  theoretical  methods  have  been  established  in  critical  comparisons  of  different 
levels  of  theory,  as  well  as  comparisons  to  experiment  for  reactions  (R4)-(R6).  In  addition, 
comparisons  between  the  BAC-N1P4  method  and  multireference  configuration  self  consistent  field 
(MCSCF)  theory  has  been  made  for  the  unimolecular  hydrogen  migration  reaction 

H0N0-»HN02  (R7) 


The  multivalent  nature  of  the  nitro  groups  makes  them  hard  to  treat  accurately  by  quantum  chemical 
methods  and  reactions  involving  nitro  groups  generally  require  a  multireference  electronic 
u  avefunction  or  a  configuration  interaction  (Cl)  calculation  for  an  adequate  description  of  the 
electronic  structure.  A  real  concern  with  the  BAC-MP4  method  has  been  that  it  is  based  upon  a 
single  reference  (IIP)  wavefunction  and  may  give  such  a  poor  zeroth  order  description  of  the 


electronic  structure  that  the  empirical  bond  additivity  correction  cannot  make  up  for  the  errors.  The 


comparisons  of  the  BAC-MP4  and  MCSCF  calculations  for  reaction  (R7)  show  excellent 
agreement  in  geometries,  frequencies,  and  qualitative  agreement  in  the  energetics.  This  give 
encouragement  for  the  use  of  the  BAC-MP4  method  for  treating  nitro  groups  which  are  very 
important  in  propellant  chemistry. 

Another  interest  that  has  developed  during  the  course  of  this  research  has  been  the 
estimation  of  kinetic  parameters  for  gas-phase  reactions.  The  ability  to  provide  reasonable 
estimates  of  reaction  rates  for  large  polyatomic  systems  is  vital  in  the  accurate  modeling  of  the 
complex  gas-phase  reaction  mechanisms.  Benson6  has  developed  methods  for  estimating  kinetic 
data  from  bond  or  group  additivity  relationships  based  on  conventional  transition  state  theory. 
These  estimates  are  obtained  from  statistical  mechanical  evaluations  of  quasithermodynamic 
parameters  based  upon  the  geometry  and  structure  of  the  saddle  point  for  the  reaction.  This 
approach  lacks  important  dynamical  effects:  (I)  the  fact  that  the  dynamical  bottleneck  controlling 
the  reaction  rate  may  be  at  a  location  different  from  the  saddle  point,  and  (2)  that  quantum 
mechanical  tunneling  makes  sizeable  contributions  to  the  thermal  rates.  These  deficiencies  have 
been  addressed  by  basing  a  thermochemical  kinetic  analysis  on  variational  transition  state  theory 
w  ith  semiclassical  adiabatic  ground-state  transmission  coefficients.  The  importance  of  these  effects 
has  been  demonstrated  and  the  development  of  a  new  database  which  incorporates  these  effects  has 
been  initiated  in  Phase  I. 

The  remainder  of  this  report  is  organized  as  follows.  Section  2  gives  a  brief  review  of  the 
theoretical  methods  used  in  this  research.  Section  3  presents  the  results  and  provides  a  discussion 
of  their  significance.  Section  4  summarizes  our  conclusions  and  outlines  the  future  directions  of 
this  research.  Two  appendices  of  papers  which  are  ready  for  submission  for  publication  are  also 
included.  Appendix  A  is  a  detailed  description  of  our  calculations  on  reactions  (R4)-(R6). 
Appendix  B  is  a  presentation  of  the  new  thermochemical  kinetic  analysis.  These  papers  include 
more  details  of  both  the  theory  and  results  for  the  interested  reader. 
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2.  Theoretical  Background 

2. 1  Variational  transition  state  theory 

Transition  state  theory  (TST)7  provides  one  of  the  most  practical  approximation  schemes 
for  calculating  equilibrium  rate  constants  for  thermally  activated  gas-phase  reactions.  The 
popularity  of  TST  over  the  last  50  years  is  largely  due  to  its  success  -  it  incorporates  the  factors 
most  crucial  in  determining  the  rates  -  and  its  simplicity  -  the  rates  are  obtained  from  equilibrium 
statistical  mechanical  calculations  and  use  only  limited  PES  information.  However,  the  utility  of 
conventional  TST  for  predicting  reaction  rates  is  limited  by  the  accuracy  of  the  approximations 
employed.  The  major  deficiencies  of  TST  are:  the  fundamental  assumption  that  if  the  reaction  gets 
to  a  critical  configuration  (the  transition  state),  then  it  will  always  proceed  to  products;  the  need  to 
incorporate  further  approximations  to  include  quantum  mechanical  tunneling  effects;  and  the 
approximate  method  bv  which  anharmonic  terms  in  the  potential  are  treated.  These  approximations 
can  lead  to  quantitatively  inaccurate  rate  constants  and  incorrect  physical  interpretations. 

Variational  transition  state  theory  (VTST)4  provides  a  convenient  framework  for 
consistently  improving  upon  the  limitations  of  conventional  TST.  The  variational  procedure  is  a 
method  to  minimize  the  errors  caused  by  the  breakdown  of  the  fundamental  assumption  of  TST. 
Conventionally  the  transition  state  is  a  dividing  surface  separating  reactants  and  products  and 
located  at  the  saddle  point  of  the  potential  energy  surface.  Significant  improvement  over 
conventional  TST  can  be  obtained  by  locating  the  dividing  surface  by  a  variational  criterion. x 
Consistent  methods  of  included  quantum  mechanical  effects  and  anharmonicity  have  been 
developed  and  extensively  tested,  and  it  has  been  found  that  the  theory  is  capable  of  giving 
quantitative  accuracy  in  the  calculation  of  thermal  rate  constants  for  atom-diatom  reactions.9' 

More  recently,  application  of  VTST  has  been  made  to  a  four-atom  system  (OH  +  H2)16  and  a 
general  computer  program  is  now  available  for  calculating  thermal  rate  constants  for  reactions 
involving  a  large  number  of  atoms. 1 7 


The  emphasis  of  this  research  has  been  the  effective  use  of  PES  information  within  the 
framework  of  VTST.  Therefore,  the  remainder  of  the  technical  discussion  is  a  brief  review  of 
some  of  the  details  of  the  theory  and  a  description  of  the  type  of  PES  information  that  is  needed  for 
the  calculations.  More  details  of  the  theory'  are  provided  in  ref.  4. 

Conventional  transition  state  theory7  reduces  the  calculation  of  the  rate  constant  to  one  of 
quasiequilibrium  statistical  mechanics:  the  equilibrium  rate  is  approximated  as  the  equilibrium  flux 
headed  towards  reactants  through  a  dividing  surface  located  at  the  saddle  point.  With  this 
approximation,  the  rate  constant  k'(T)  for  temperature  T  takes  on  the  simple  textbook  form 


k  ( T) 


kBT 

h 


Q'(T) 

Ok(T) 


(1) 


where  kQ  is  Boltzmann's  constant,  h  is  Planck's  constant,  Q'(T)  is  the  partition  function  for  the 
bound  degrees  of  freedom  at  the  saddle  point,  O^T)  is  the  reactants  partition  function,  and  Vf  is 
the  value  of  the  potential  at  the  saddle  point.  Thus  conventional  transition  state  theory-  requires 
information  about  the  potential  energy  surf  ace  only  in  the  saddle  point  and  reactant  regions.  If  the 
partition  functions  are  computed  using  a  harmonic  approximation  then  the  matrix  of  second 
derivatives  (the  Hessian  matrix i  suffices. 

In  variational  transition  state  theory  the  dividing  surface  is  viewed  as  a  tentative  dynamical 
bottleneck  to  flux,  and  the  best  bottleneck  (the  dividing  surface  allowing  the  least  flow  of  flux)  is 
located  by  a  variational  procedure.  The  generalized  expression  for  the  thermal  rate  constant  is 
given  as  :>  function  of  the  location  s  of  the  dynamical  bottleneck  along  the  reaction  coordinate 


k  T  ,  XiT  ...  , 

C)T  K  B  Q  (  T.  s ) 

k  (T.  s)  =  —  exP(  V  Ml.p(s)  /  k  BT) 


where  (/’  r(T.s>  is  the  generalized  partition  function  for  the  bound  degrees  of  freedom  orthogonal 
to  the  reaction  path  at  s,  and  VN1[:[>(s)  is  the  value  of  the  potential  along  the  reaction  path  at  s.  One 


version  of  variational  transition  state  theory,  the  canonical  variational  theory  (CVT),  results  from 
minimizing  eq.  (2)  with  respect  to  s 


.  CVT  .  .  GT 

k  (T)  =  min  k  (T,  s) 

s 


(3) 


The  improved  canonical  variational  theory  (ICVT)1"  also  variationally  optimizes  the  location  of  the 
transition  state  dividing  surface  for  a  given  temperature,  but  provides  an  improved  treatment  of 
threshold  energies  by  using  an  ensemble  which  removes  energies  below  the  ground-state  adiabatic 
threshold.  To  compute  the  rate  constant  using  either  the  canonical  or  improved  canonical 
variational  theory  more  information  about  the  potential  energy  surface  is  required  than  fora 
conventional  transition  state  theory'  calculation;  information  about  the  potential  in  a  region  around 
the  reaction  path  is  required. 

For  most  chemical  reactions  of  interest,  especially  those  involving  hydrogen  atoms,  it  is 
imperative  to  include  quantum  mechanical  effects  in  the  theory.  Quantum  mechanical  effects  are 
included  in  VTST  in  an  ad  hoc  manner:  bound  modes  are  quantized  by  using  quantum  mechanical 
partition  functions  in  eq.  (2)  and  quantum  mechanical  effects  on  the  reaction  coordinate  motion  are 
included  by  a  multiplicative  factor,  the  transmission  coefficient  kx(T).  The  quantized  rate  constant 
is  given  by 


ICVT  /  X 


<  T)  =  KX(  T) 


ICVT 

k  (T) 


(4) 


where  the  superscript  X  specifies  the  yet  undefined  method  used  in  calculating  the  transmission 
coefficient.  The  transmission  factor  includes  the  effect  of  nonclassical  reflection  for  energies  above 
the  classical  barrier  as  well  as  the  effect  of  quantum  mechanical  tunneling.  The  Wigner  correction 
factor1*  has  been  verv  popular  for  estimating  quantum  mechanical  tunneling  effects  because  it 
requires  knowledge  of  the  shape  of  the  classical  barrier  near  the  saddle  point;  this  type  of 
information  can  be  obtained  from  the  negative  eigenvalue  of  the  Hessian  at  the  saddle  point 
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Unfortunately,  the  approximations  of  this  method  are  rarely  valid  for  chemical  reactions,  i.e., 
w  hen  tunneling  is  important  it  occurs  over  distances  much  longer  than  that  represented  by  a 
quadratic  representation  of  the  potential  near  the  saddle  point,  and  a  more  global  description  of 
tunneling  is  required. 

More  successful  transmission  coefficients  are  based  upon  the  vibrationally  and  rotationally 
adiabatic  approximation.19’-0  In  this  approximation  the  difficult  multidimensional  scattering 
problem  is  transformed  into  the  much  easier  problem  of  scattering  on  a  one-dimensional  effective 
potential.  The  vibrationally-rotationallv  adiabatic  potential  is  given  by 

vaU«)  =  vMEP(s)  +  eST(s)  (5) 

where  a  is  a  collective  index  of  the  quantum  numbers  for  the  bound  modes  and  e^T(s)  is  the 
hound  energy  level  for  state  a  at  the  generalized  transition  state  located  at  s  along  the  reaction  path. 
The  adiabatic  approximation  is  made  in  a  curvilinear  coordinate  system  and  although  the  potential 
term  is  very  simple,  the  kinetic  energy  term  is  complicated  by  factors  dependent  upon  the  curvature 
of  the  reaction  path  (the  MEP).  The  method  of  treating  the  reaction-path-curvature  terms  is 
ambiguous  and  several  approaches  have  been  tried.  The  best  methods  treat  tunneling  along  a  path 
which  depends  upon  the  reaction  path  curvature.  An  example  is  the  small-curvature  semiclassical 
ground  state  (SCSAG)  tunneling  method  which  is  based  upon  the  ground-state  adiabatic  potential 
and  in  which  the  tunneling  probabilities  are  computed  by  the  WKB  method.20  This  method  works 
well  for  systems  in  which  the  reaction  path  curvature  is  not  too  large.  Note  that  for  this  tunneling 
method  the  type  of  PES  information  needed  is  nearly  identical  to  that  needed  for  the  variational 
calculation:  the  only  additional  piece  of  information  necessary  is  the  curvature  of  the  reaction  path 
w  Inch  can  be  obtained  from  the  first  and  second  derivatives  of  the  potential  along  the  reaction  path. 

In  the  present  work  the  ICVT/SCSAG  method  is  applied  using  limited  information  about 
potential  energy  surface:  the  information  necessary  to  carry  out  these  calculations  is  the 
■  as  :;.t:  along  the  minimum  energy  path  (MEP:  the  minimum  energy  path  is  the  path  of  steepest 
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docents  in  a  mass  scaled  coordinate  system  from  the  saddle  point  to  both  reactants  and  products), 
the  matrix  of  force  constants  along  the  minimum  energy  path,  and  the  curvature  of  the  reaction 
path.  The  method  directly  uses  the  results  of  quantum  chemistry  calculations  of  information  about 
the  potential  energy  in  the  form  of  the  energy  and  its  derivatives.  The  development  of  quantum 
chemistry  methods  to  analytically  evaluate  the  derivatives21  has  greatly  reduced  the  computational 
effort  and  such  methods  are  used  when  available.  Typically  the  potential  information  is  obtained 
on  a  sparse  grid  along  the  minimum  energy  path  and  methods  have  been  developed  to  interpolate 
the  necessary  parameters  of  the  reaction-path  Hamiltonian  (RPH).1 


2.2  Quantum  chemistry 

Today,  most  quantitative  theoretical  studies  of  intermolecular  forces  and  their  related 
potential  energy  surfaces  are  based  on  ab  initio  methods,  which  can  be  pursued  to  increasing 
levels  of  accuracy,  thereby  providing  useful  estimates  of  the  errors  at  preceding  levels. 

Historically,  the  majority  of  calculations  have  determined  the  ground  electronic  states  of  atomic  and 
molecular  systems.  The  most  well-known  ab  initio  method  is  the  LCAO-SCF  or  Hartree-Fock 
(HF)  method.2"1  In  this  procedure  the  wavefunction  for  the  system  is  expressed  as  a  determinant 
of  one-electron  orbitals  formed  from  linear  combinations  of  orbitals  centered  on  each  of  the  nuclei. 
The  Hartree-Fock  solution  is  a  "variational"  one:  i.e.,  for  a  given  set  of  atomic  orbitals  ("basis 
set"),  it  yields  the  lowest  energy  possible  for  a  single-determinant  wavefunction.  In  this  type  of 
calculation  there  can  be  two  sources  of  error:  that  from  the  use  of  a  finite  basis  set  and  that  from  the 
neglect  of  correlation.  Experience  has  show-n  that  the  accurate  determination  of  properties, 
particularly  bond  strengths,  even  for  singly  bonded  substituents,  requires  that  both  types  of  errors 
lx*  eliminated  to  the  maximum  extent  possible. 

Theoretically,  the  correlation  energy  is  defined  to  be  the  difference  between  the  Hartree- 
Fock  energy  and  the  exact  ( nonrelati vistic )  energy  within  the  basis.24  When  a  method  corrects  for 
the  approximations  in  the  independent  particle  model  by  including  components  from  other  electron 
configurations  in  the  wavefunction,  correlation  is  accounted  for  and  the  total  electronic  energy 
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obtained  will  be  lower  than  the  HF  energy.  The  conventional  method  of  obtaining  a  correlated 
wavefunction  is  the  configuration  interaction  (Cl)  method.25  When  increasing  numbers  of 
determinants  of  orbitals  representing  states  excited  relative  to  the  SCF  determinant  are  included  in 
the  description  of  the  wavefunction,  the  Cl  method  converges  monotonically  (from  above)  to  the 
basis-set-limited  exact  electronic  energy  (i.e.,  the  "full  Cl"  limit).  The  convergence  with  respect  to 
the  basis  set  can  be  assessed  by  examining  the  changes  in  calculated  properties  as  basis  functions 
are  added.  It  is  particularly  important  that  the  errors  in  different  regions  of  a  potential  energy 
surface  remain  fairly  constant:  otherwise,  spurious  barriers  or  other  features  may  appear. 

Electron  correlation  may  be  classified  as  dynamical  or  nondynamical.26  Nondynamical 
correlation  errors  are  large  for  systems  which  are  inadequately  described  by  a  single  determinant. 
This  case  arises  when  two  or  more  configurations  are  nearly  degenerate  in  energy  or  when 
alternative  arrangements  of  electrons  in  partially  filled  subshells  are  required  for  proper 
dissociation.  One  method  particularly  appropriate  for  obtaining  nondynamical  correlation  is  the 
multiconfiguration-SCF  (MCSCF)  method.27  In  this  procedure,  both  the  orbital  coefficients  and 
the  correlation  coefficients  within  a  limited  Cl  space  are  varied  and  iterated  to  convergence.  For  a 
given  configuration  space,  the  MCSCF  orbitals  are  the  optimal  ones.  Dynamical  correlation  is 
obtained  in  this  procedure  by  a  Cl  using  the  optimized  orbitals  to  give  the  MCSCF-CI28  method. 

An  alternate  to  the  Cl  method  is  based  upon  perturbation  theory  such  as  diagrammatic 
many-body  perturbation  theory  (MBPT)29  or  Moller-Plesset  perturbation  theory.30  As  higher- 
order  perturbations  are  included,  the  energy  (for  systems  for  which  a  single  determinantal  starting 
point  is  adequate)  converges  to  the  full  Cl  limit,  although  not  usually  monotonically  since  the 
MBPT  energy  does  not  include  certain  spurious  positive  terms  present  in  nonfull  Cl  energies  that 
are  cancelled  as  the  full  Cl  is  reached  and  are  responsible  for  the  monotonic  convergence. 
However,  recognition  of  these  terms  makes  it  possible  to  approximately  correct  and  improve  the 
nonfull  Cl  energies.31  For  ground  state  thermochemistry  and  properties  at  equilibrium  geometries 
and  transition  states,  the  perturbation  based  codes  (e.g.,  GAUSSIAN  8232)  offer  superior 
performance  over  MCSCF-CI  calculations  in  many  instances. 


In  the  present  work  we  attempt  to  keep  the  computational  effort  affordable  by  using 
moderate  sized  basis  sets  and  partial  descriptions  of  the  electron  correlation.  These  give  results  of 
modest  accuracy  which  are  improved  by  empirical  corrections  using  the  bond  additivity  correction 
(BAC)  method.-''  This  approach  begins  with  a  Hartree-Fock  (HF)  optimization  of  critical 
geometries  (such  as  equilibrium  geometries  and  transition  states),  followed  by  higher  quality 
calculations  of  the  energy  using  full  fourth-order  Moller-Plesset  perturbation  theory-  (MP4),  and 
finally  applying  an  empirical  bond  additivity  correction.  This  method  is  denoted  HF/BAC-MP4  to 
distinguish  it  from  other  possibilities  discussed  below.  The  correction  is  given  in  terms  of  an 
adjustment  to  the  total  energy  for  each  bond  in  the  molecule.  The  bond  correction  dies  off 
exponentially  with  the  bond  length.  The  parameters  of  the  BAC-MP4  method  are  fitted  to  a 
database  of  experimental  heats  of  formation  of  atomic  and  molecular  species  in  a  least  squares 
sense.  These  parameters  will  change  with  the  level  of  the  ab  initio  calculation  and  in  the  limit  of  a 
complete  basis  set  and  full  Cl  the  correction  goes  to  zero.  The  strengths  of  the  BAC-MP4  method 
are:  ( 1 )  it  is  based  upon  ab  initio  methods  and  therefore  incorporates  a  qualitatively  correct  picture 
of  the  multicenter  nature  of  the  molecular  bonding,  (2)  it  can  be  systematically  improved  by 
improving  the  level  of  the  underlying  ab  initio  calculation,  and  (3)  the  corrections  are  estimated 
from  a  large  database  of  experimental  heats  of  formation.  For  corrections  to  energies  at 
equilibrium  geometries  the  BAC  method  can  be  viewed  as  an  interpolation  technique,  whereas  for 
transition  states  and  other  geometries  in  the  strong  interaction  region  for  which  experimental  data  is 
not  available  the  method  is  more  an  extrapolation  technique  and  must  be  more  critically  tested. 

A  major  concern  in  using  the  HF/BAC-MP4  method  is  the  accuracy  of  the  computed 
potential  energy  surface  information.  One  test  is  to  improve  the  quality  in  the  underlying  ab  initio 
calculation  to  decrease  the  magnitude  of  the  empirical  bond  additivity  correction  and  test  the  effect 
on  the  computed  rate  constant.  As  a  first  step  in  this  direction,  we  have  reevaluated  the  potential 
information  using  the  MP2/BAC-MP4  method.  In  this  method  the  geometry  of  the  equilibrium 
structure  or  transition  state  is  determined  using  a  full  second-order  Moller-Plesset  perturbation 
theory  (MP2)  calculation.  The  matrix  of  second  derivatives  is  also  calculated  at  the  MP2  level  at 


the  same  geometry.  The  minimum  energy  path  is  located  by  following  the  path  of  steepest 
descents  using  the  gradient  vectors  and  Hessian  matrices  computed  at  the  MP2  level.  For  each 
point  along  the  minimum  energy,  the  energy  is  then  computed  by  the  MP4  method  followed  by  a 
bond  additive  correction.  Because  the  geometries  and  vibrational  frequencies  are  different  at  the 
I  IF  and  MP2  levels,  the  BAC  correction  is  also  different  and  the  parameters  are  refitted  to 
experimental  heats  of  formation  as  was  done  for  the  HF/BAC-MP4  method. 

One  further  step  that  has  been  considered  is  the  MP4/BAC-MP4  method  in  which  the 
energy  and  first  and  second  derivatives  along  the  entire  minimum  energy  path  are  calculated  by  the 
MP4  method.  The  bond  additivity  correction  is  then  made  to  the  MP4  energy  at  each  point  along 
the  minimum  energy  path.  This  method  has  not  been  fully  implemented,  but  by  observing  the 
change  of  the  geometries  and  frequencies  in  going  from  the  MP2/BAC-MP4  method  to  the 
MP4/BAC-MP4  method,  we  obtain  an  estimate  of  the  sensitivity  of  the  calculations  to  even  higher 
levels  of  the  underlying  ab  initio  electronic  structure  calculations.  Similar  comparisons  have  been 
made  between  the  HF/BAC-MP4  method  and  MCSCF  calculations  for  which  a  bond  additivity 
correction  is  not  available.  These  comparisons  give  an  estimate  of  the  importance  of  using  a 
multireference  wavefunction  for  optimizing  geometries  and  obtaining  frequencies,  and  provide  a 
good  measure  of  the  accuracy  of  the  HF-based  pertubative  approaches. 

2.3  Thermochemical  kinetic  analysis 

Our  ability  to  address  reactions  involving  large  polyatomic  molecules  such  as  found  in  the 
initial  stages  of  nitrocubane  decomposition  will  rely  to  a  great  extent  on  developing  methods  to 
estimate  Arrhenius  parameters  from  bond  additivity  or  group  additivity  relationships.  Statistical 
mechanics  provides  the  prescription  for  using  information  about  the  structure  and  vibrational 
frequencies  of  bound  chemical  species  to  directly  calculate  thermochemical  data  such  as  heats  of 
formation  of  stable  compounds  and  heat  release  and  equilibrium  constants  for  chemical  reactions.  ^ 
The  structural  and  energetic  information  can  be  reduced  into  thermodynamic  quantities,  the 


enthalpy  and  entropy,  which  can  be  accurately  approximated  using  bond  additivity  and  group 
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additivity  relationships.  Using  the  quasithermodynamic  formulation  of  conventional  transition 
state  theory6,7  these  ideas  have  been  extended  to  the  prediction  of  reaction  rate  constants  by 
estimating  the  enthalpy  and  entropy  of  the  activated  complex  from  knowledge  of  the  structure  and 
s  ibrational  frequencies  at  the  saddle  point  of  the  reaction.  This  approach  neglects  important 
contributions  to  the  thermal  rate  from  variational  optimization  of  the  transition  state  and  quantum 
mechanical  tunneling.  Recently,  a  new  partitioning  of  the  thermal  rate  constant  has  been  produced 
v.  hich  allow  s  analysis  of  the  rates  in  terms  of  quasithermodynamic  quantities  which  include  these 
important  dynamical  effects.'4 

In  traditional  thermochemical  kinetic  analysis  the  enthalpy  and  entropy  of  activation  are 
estimated  from  properties  at  a  single  location  along  the  minimum  energy  path  -  the  saddle  point. 
Basing  the  thermochemical  kinetic  analysis  upon  variational  transition  state  theory  presents  a 
problem  because  the  role  of  the  saddle  point  is  replaced  by  that  of  the  dynamical  bottleneck  w  hich 

moves  with  temperature.  At  0  K  it  occurs  at  the  maximum  of  the  ground-state  adiabatic  potential 
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curve  s'  „ 1  and  past  experience  has  shown  that  the  dynamical  bottleneck  for  temperatures  up  to 

about  6(H)- 1000  K  tend  to  be  in  the  vicinity  of  this  adiabatic  maximum.  Therefore,  it  is  more 

reasonable  to  base  the  thermochemical  kinetic  analysis  upon  the  generalized  transition  state  rate 

AG 

constant  given  by  eq.  (2  )  evaluated  at  s  ,  instead  of  the  saddle  point  location  of  s=0.  This  takes 

into  account  most  of  the  variational  effects  but  still  neglects  thermal  fluctuations  of  the  dynamical 
AG 

bottleneck  from  s'  „  which  become  more  important  at  higher  temperatures,  and  also  neglects 

quantum  mechanical  tunneling  effects  w  hich  are  important  at  low  temperatures.  These  effects  are 
also  included  by  partitioning  the  rate  constant  |eq.(4)|  as  follows 


,  1CVT/X  ... 
k  (  I  )  = 


(T)  k°r(T,  s')G> 


where  K1)k,m)(T>  accounts  for  finite-temperature  deviations  due  to  the  temperature  dependence  of 
the  improved  canonical  variational  transition  state  and  is  given  by  the  expression 
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k'CVT|T>  =  NlK-n„('D  kflT(T.SA.°) 

Ihe  generalized  transition  state  theor>’  rate  constant  k^T(T,sA^)  is  obtained  from  the  properties  of  a 

single  generalized  transition  state,  "a  substance",  and  is  factored  in  the  usual  manner  in  terms  of 
substantial  activation  parameters,  AsS<,j.  and  ASH^, 


kG1(T,  s)  = 


YtcV"  exP(AsST/  R)cxp(- ASH°/RT) 
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v*.  here  O'  is  the  concentration  in  the  standard  state  and  An*  is  the  stoichiometric  change  in  the 
number  of  moles  in  passing  from  the  reactants  to  the  transition  state  (An*=-1  for  reactions  studied 
here).  For  an  ide; '  gas  at  a  standard  state  of  P°  =  1  atm,  the  concentration  of  the  standard  state  is 
given  by  P°/RT,  and  the  substantial  enthalpy  and  entropy  of  activation  can  be  expressed  in  terms 
of  the  reaction  rate  constant  by 


A  H1’  =  RT  2  ^-1^-  -  2RT 
s  I  d  1 
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where  \A  is  Avogadro's  number,  and  u  is  one  for  k(T)  in  molar  units  and  zero  for  k(T)  in 
molecular  units. 

The  transmission  coefficients  K^rmCD  and  )  depend  upon  more  extended 

regions  of  the  potential  energy  surface  rather  than  a  single  transition  state  location  but  we  still  factor 
them  in  a  manner  analogous  to  eq.  (K),  in  terms  of  nonsubstantial  activation  parameters,  AnS  y 
and  AnIl". 

M  T>  =  exp  (  A  nSy  /  R)exp(-  A„Il"  /RT)  (11) 


16 


V  •/  A  V  ■/  V 

i  -vfc  A  A  JW  -  ^  .  *- 


■VA  .  /-VA  A  A  A  A  Ay/i  /.  A  A  /^\  A  .V  A  A  ^  AA  AA  A.  AjlA 


(12) 


t 


s 


I 

► 

I 

t 

t 

0 


I 

* 

» 

i 


* 

IS 

S, 

K. 

K. 

S. 

S. 

K. 

I 


w  here 


mT)=  kSCSaG(T)  k  ,  (T) 

tiie  mi 


and  the  nonsubstantial  enthalpy  and  entropy  of  activation  are  obtained  from 


=  rt  2 


AnS°  =  RT  d["{K  +  R  In  k(T) 


(13) 


(14) 


Since  AsS°  and  AJH^  are  quasithermodynamic  parameters,  their  temperature  dependences  can 
be  obtained  from  a  single  temperature-dependent  substantial  heat  capacity  of  activation  A  SC ° 


Vt-Vt  +it  (ASCP/T)  dT 


and 


AsH”=Asn"  +^AsCpdT 

1  ()  o 


(18) 


(19) 


The  parameters  AnS  and  A n H  T  are  purely  phenomenological;  however,  because  of  their 

definition  in  eqs.  ( 13)  and  ( 14),  their  temperature  dependences  can  also  be  obtained  from  a  single 
heat  capacity  analog  AnC ;  given  by  expressions  analogous  to  eqs.  (15)  and  (16). 

This  new  thermochemical  kinetic  analysis  separates  contributions  from  a  single  transition 
state  structure  from  contributions  due  to  thermal  and  quantum  mechanical  tunneling  which  depend 
on  more  global  properties  of  the  potential  energy  surface.  The  former  can  be  estimated  using  bond 
additivity  or  group  additivity  relationships,  but  the  new  schemes  are  needed  for  estimating  the 
latter. 
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3.  Results  and  Discussion 


3. 1  Validation  of  theoretical  methods 

3.1.1  Application  to  II  +  NH,  ^  Ht  +  NH2  and  D  +  ND3  ->  D,  +  ND, 

Greater  detail  and  discussion  of  the  results  for  these  reactions  as  well  as  a  prediction  of  the 

|  rate  constants  for  the  reverse  of  reaction  (R6)  is  presented  in  appendix  A.  The  results  are 

summarized  in  this  section  with  an  emphasis  on  the  validation  of  the  HF/BAC-MP4  method  and 
the  computational  effort  required  for  these  calculations. 

j  First  the  saddle  point  geometries  are  compared  at  the  HF,  MP2,  and  MP4  levels  for  the 

reaction  of  H  w  ith  NH3.  At  the  saddle  point  the  geometry  is 


where  the  bond  lengths  and  angles  are  given  in  table  1. 


Table  1.  Saddle  point  geometry  for  the  H  +  NH3  — >  FU 

+  NH2  reaction. 

HF 

MP2 

MP4 

rnh/A> 

1.244 

1.305 

1.277 

RN[I:(A) 

1.010 

1.023 

1.026 

Rnh/A) 

1.010 

1.023 

1.026 

RH1I!./'°X) 

0.949 

0.869 

0.893 

^'HjH /degrees) 

100.5 

99.1 

99.3 

100.5 

99.1 

99.3 

NH  Agrees) 

105.4 

103.5 

103.3 

V. 


I 


y. 


i 

v 

V. 

s'. 

V 

V 

s\ 


geometries.  The  largest  changes  are  for  the  NH,  and  HjH4  bond  lengths  which  are  most  actively 

involved  in  the  rearrangement.  Note  that  the  geometries  do  not  change  monotonically  with 
improvement  in  the  level  of  treatment  but  the  HF  values  for  RN1(i  and  tire  closer  to  the 

MP4  values  than  are  the  MP2  values.  The  NH]  and  H]H4  bond  distances  computed  at  the  HF  and 

MP2  levels  are  show  n  m  figure  1  tor  points  along  the  minimum  energy  path.  The  agreement 

between  the  two  methods  for  the  whole  minimum  energy  path  is  as  good  as  or  better  than  the 

agreement  seen  at  the  saddle  point  . 

The  saddle  point  frequencies  tire  compared  at  the  HF.  MP2.  and  MP4  levels  for  the  reaction 
of  H  w  ith  N1 1-  in  table  2.  The  imaginary  frequencies  at  the  top  of  each  column  correspond  to  the 

unbound  mode  for  motion  along  the  reaction  coordinate.  Note  that  of  all  the  frequencies  this 
imaginary  one  changes  the  most  as  the  level  of  theory  is  improved.  The  harmonic  zero  point 
energy  given  at  the  bottom  of  each  column  is  one-half  the  sum  of  the  frequencies  for  all  bound 
modes.  Although  the  frequencies  for  some  individual  modes  do  change  by  as  much  as  20%,  the 
changes  in  zero  point  energy  are  much  smaller.  A  change  of  1 10  cm'1  in  the  zero  point  energy  (a 
typical  difference  between  the  HF  and  MP2  values)  will  shift  the  adiabatic  potential  curve  by  only 
0.3  kcal/mol. 
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Figure  1.  NH  and  HH  bond  distances  along  the  minimum  energy  path  for  the  reaction  H  + 

— >  H-,  +  NUi.at  two  levels  of  ah  initio  electronic  structure  theory  -  Hartree-Fock  (*)  and  second 
order  Moller  Plesset  perturbation  theory  (•) 


C Vvs- .'■j.'v'v  v'v  ’  '  V  V. 
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Table  2.  Saddle  point  harmonic  frequencies  zero  point  energy  (cm'1)  for  the  H  +  NH-,  — >  H-,  + 
N1 1 -..reaction. 


Zero  point  enercv 


HF 

MP2 

MP4 

283  vi 

202  li 

1905i 

723 

718 

726 

732 

730 

901 

1258 

1 146 

1186 

1485 

1315 

1347 

1541 

1610 

1629 

1731 

2077 

1921 

3647 

3489 

3497 

3750 

3606 

3553 

7434 

7345 

7381 

The  bond  additivity  correction  is  based  upon  information  about  the  potential  energy  surface 
only  through  quadratic  terms  (a  harmonic  treatment)  but  is  designed  to  fit  the  experimental  heats  of 
formation  of  bound  molecules  which  include  the  effects  of  anharmonicity.  Therefore,  the  BAC 
method  has  anharmonic  effects  built  into  it  in  an  averaged  sense,  not  on  a  mode-by-mode  basis. 

The  harmonic  frequencies  computed  at  the  HP  level  are  known  to  be  systematically  high  and  the 
BAC  method  corrects  for  this  discrepancy  as  well  as  including  the  effects  of  anharmonicity.  The 
relatively  small  differences  in  the  zero  point  energy  w  ith  improved  level  of  theory  is  encouraging 
since  it  means  that  smaller  changes  will  have  to  be  made  in  the  bond  additivity  correction  to  account 
for  these  changes  in  frequencies. 

A  comparison  of  the  potential  energy  curve  V\1f.p(s)  and  the  ground-state  adiabatic 
potential  curves  VA(s)  for  reaction  (R4)  using  the  HF/BAC-MP4  and  MP2/BAC-MP4  methods  is 


cthod  umik  the  potential  information  obtained  by  the  IIF/BAC-MP4  and  MP2/BAC-MP4 


methixls  is  given  in  table  3  and  figure  3. 


Table  3.  Comparison  of  ICVT/SCSAG  reaction  rates  (in  units  of  cnv1molecule'1s'1 )  computed 
using  two  different  sets  of  information  about  the  potential  energy  surface. _ 


T.K 

HF/BAC-MP4 

MP2/BAC-MP4 

200 

9.9x10- 24 

3.2x1  O'24 

300 

1.8x10'“° 

1 .2x1 0"20 

400 

2.0x1  O' 1S 

1.6xl0‘1K 

600 

4.5x10- 16 

4.0.x  10' 16 

1000 

6.2x1  O' 14 

6.2x  Kb 14 

1500 

l.lxlO" 12 

l.lxlO"12 

2400 

6.9x10- 12 

8. 1  x  1 0" 1 2 

Also  included  in  figure  3  is  a  comparison  with  the  experimental  results  of  Marshall  and  Fontijn.-0 
The  potential  along  the  minimum  energy  path  and  the  adiabatic  potential  curves  show  very'  similar 
shapes  for  the  two  methods,  although  the  MP2/BAC-MP4  results  are  shifted  slightly  to  the  left. 
Tire  classical  barrier  heights  are  very  similar.  17.1  keal/mol  for  HF/BAC-MP4  compared  to  16.7 
keal/mol  for  MP2/BAC-MP4.  and  the  adiabatic  barriers  are  15.5  and  15.7  keal/mol  above  the 
reactant  zero  point  energies  for  the  IIF/BAC-MP4  and  MP2BAC-MP4  methods,  respectively. 

The  rate  constants  computed  using  the  two  sets  of  potential  information  agree  well  for 
temperatures  .ilvve  300  K.  At  300  K  and  below  tunneling  is  very  important  and  the  rates  tire 
sensitive  to  finer  details  of  the  adiabatic  potential  curves.  At  300  K.  the  rates  computed  using  the 
two  sets  nt  potential  intomiation  differ  by  only  50C  but  at  200  K.  the  difference  increases  to  a 
-,u  tor  ot  :  ! .  fhe  overall  agreement  with  the  experimental  results  is  also  very  good  except  at  the 


Figure  2.  Comparison  of  potential  energy  along  the  minimum  energy  path  I  VMFf>(s)|  and  ground- 

state  adiabatic  potential  |VA(s)|  as  a  function  of  reaction  coordinate  s  computed  using  the  HF/BAC 
MP4  method  (solid  curves)  and  the  MP2/BAC-MP4  method  (dashed  curves)  for  the  reaction  1 1  + 
Ml,  —>11-,  +  Ml-,. 


WT 


v\ 
V*  ‘ 


O  experiment 
-  HF/BAC-MP4 
MP2/BAC-MP4 


1000/T  (IT1 ) 


I  mure  V  Rato  constants  as  a  function  of  temperature  for  the  reaction  II  +  NH3  —»  +  N1U 

computed  using  the  improved  canonical  variational  theory  with  small-curvaturc  scmiclassical 
.id;.il\ttic  ground -state  transmission  coefficient  from  potential  information  obtained  using  the 
I II  PA  (MP4  method  (solid  curves)  and  the  MP2/BAC-MP4  method  (dashed  curves). 

(  .>mpari \on  is  also  made  with  the  experimental  rate  constants  of  P.  Marshall  and  A.  Fontijn.  ref 

1  '5 1.  denoted  Mb  (  •  ). 


lowest  temperatures  which  is  discussed  in  more  detail  below.  Comparison  can  also  be  made  to 
least-squares  Arrhenius  fits  to  experimental  results  for  the  temperature  ranges  from  500  to  1 140 
K.'-'  90S  to  1777  K,-6  and  673  to  1003  Kf7  These  fits  give  activation  energies  of  17.2±0.8, 
!6.0±0.2.  and  14.6±1.0  kcal/mol,  respectively.  Activation  energies  derived  from  the  computed 
rates  show  an  interesting  behavior  of  initially  increasing  w'ith  increasing  temperature  but  turning 
around  and  decreasing  at  even  higher  temperature;  e.g.,  the  computed  activation  energies  are  12.9, 
14.7.  17.1.  and  14.7  kcal/mol  over  the  temperature  ranges  400-600  K,  600-1000  K,  1000-1500 
K.  and  1500-2400  K  respectively.  In  general  the  agreement  between  the  experimental  and 
computed  activation  energies  is  excellent:  however,  the  theoretical  results  systematically 
underestimate  the  experimental  reaction  rate  indicating  a  systematic  overestimate  in  the  free  energy 
of  activation.  The  underestimation  of  the  rate  may  be  attributed  to  a  barrier  height  which  is  too 
large:  however,  lowering  just  the  barrier  will  tend  to  decrease  the  activation  energy  which  is 
inconsistent  w  ith  the  experiments.  Another  possibility  is  that  the  energetics  are  approximately 
correct  (i.e..  the  ground-state  adiabatic  barrier  is  about  right)  and  that  the  entropic  factors  controlled 
by  the  higher  lying  vibrational  levels  are  not  accurately  described. 

Reaction  rates  computed  by  the  the  ICVT/SCSAG  method  with  the  HF/MP4-BAC  potential 
information  for  reaction  (R5)  are  compared  with  the  experimental  ones  of  Hack  et.  al ,77  and 
Demissy  and  Lesdaux-1**  in  figure  4.  The  agreement  w  ith  the  low'  temperature  results  of  Demissy 
and  I.esclaux  is  excellent  (differences  of  18-40%)  but  at  higher  temperatures  the  computed  rates  are 
factors  of  2.2  to  2.7  lower  than  the  results  of  Hack  et.  al.  The  discrepancy  at  the  higher 
temperatures  can  be  attributed  to  uncertainties  in  the  heat  of  formation  of  NH-,.  At  low-er 
temperatures  the  computed  heat  of  formation  of  NH?  agrees  well  w  ith  the  experimental  values  and 
we  expect  less  uncertainty  in  the  computed  rate  constant. 

Reaction  rates  computed  by  the  ICVT/SCSAG  method  using  the  HF/RAC-MP4  potential 
information  for  reaction  (R6)  are  compared  with  experimental  one  of  Marshall  and  Fontijn™  in 
tie  are  5.  Once  again  the  agreement  is  excellent,  although  the  theoretical  results  systematically 
underestimate  the  reaction  rates.  This  underestimate  is  consistent  with  our  observations  for  the  H 


molecule 


ICVT/SCSAG 


1000/T  (IT1 ) 


Figure  4.  Comparison  of  computed  and  experimental  rate  constants  for  the  reaction  +  NFF,  — 
1 1  +  N if v  The  experimental  results  are  the  Arrhenius  fit  of  W.  Hack,  P.  Rouveirolles,  and  H. 
Gg.  Wagner,  ref.  (36),  denoted  I IRW  (short  dashed  curve)  and  the  Arrhenius  fit  of  M.  Demissy 
and  R.  Lesclaux,  ref.  (38),  denoted  DL  (long  dashed  curve).  The  solid  line  is  the  result  of  the 
improved  canonical  variational  theory  with  small-curvature  semiclassical  adiabatic  ground-state 
transmission  coefficient  (ICVT/SCSAG)  from  potential  information  obtained  using  the  HF/BAC- 
MP4  method. 
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ICVT/SCSAG 
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Figure  5.  Comparison  of  computed  and  experimental  rate  constants  for  the  reaction  D  +  ND^  — » 
I)-,  +  ND-,.  The  experimental  results  are  from  P.  Marshall  and  A.  Fontijn,  ref.  (39),  denoted  MF 

( • ).  The  solid  line  is  the  result  of  the  improved  canonical  variational  theory  with  small-curvature 
semiclassical  adiabatic  ground-state  transmission  coefficient  (ICVT/SCSAG)  from  potential 
information  obtained  using  the  HF/BAC-MP4  method. 


+  Nil,  — >  1 U  +  NH-,  reaction  that  some  feature  of  the  potential  energy  surface  is  giving  rise  to  a 
free  energy  of  activation  which  is  too  low  .  Compared  to  reaction  (41),  tunneling  is  less  important 
for  the  deuterated  reaction.  For  example,  at  300  K,  the  1CVT/SCSAG  transmission  coefficient  is 
44.5  for  reaction  (R4)  and  only  8.5  for  (R5),  and  at  600  K  these  decrease  to  2.4  and  1.6, 
respectively.  As  for  reaction  (R4),  the  good  agreement  between  the  computed  and  experimental 
actis  ation  energy  indicates  that  the  reaction  energetics  are  approximately  correct  but  the  entropic 
factors  controlled  by  higher  lying  vibrational  levels  leads  to  the  underestimate  of  the  rates. 

For  the  three  reactions  for  which  experimental  results  are  available,  the  agreement  between 
theory  and  experiment  is  very'  good  except  in  comparison  with  the  experimental  results  of  Marshall 
and  Fontijn-5-5  for  reaction  (R4)  at  temperatures  below  660  K.  Below  this  temperature  the 
experimental  results  show  much  more  curvature  in  the  Arrhenius  plot  than  the  theoretical  rates  and 
at  5(X)  K  the  experimental  rates  are  higher  than  the  theoretical  ones  by  about  a  factor  of  4.  Such 
curvature  in  the  Arrhenius  plot  is  often  attributed  to  quantum  mechanical  tunneling,  but  in  the 
present  calculations  the  large  amount  of  tunneling  necessary  to  reproduce  the  experimental 
curvature  is  inconsistent  with  the  excellent  agreement  between  theory  and  experiment  for 
temperatures  above  660  K.  For  example,  at  660  K  the  SCSAG  tunneling  correction  factor  is  about 
2  and  gives  good  agreement  with  experiment.  At  500  K  the  SCSAG  tunneling  correction  factor  is 
3.6  but  would  have  to  be  about  15  to  agree  with  experiment.  To  reproduce  this  sharp  increase  in 
tunneling  over  such  a  small  change  in  temperature  w-ould  require  a  radically  different  shape  of  the 
adiabatic  potential  which  is  inconsistent  with  the  reliability  established  in  previous  tests  and  uses  of 
the  methods. 

Another  indication  of  the  inconsistency  of  these  low  temperature  results  is  given  by  the  low 
temperature  results  for  reaction  (R5)  w  hich  are  related  to  those  of  reaction  (R4)  by  detailed  balance. 
Knowledge  of  the  equilibrium  constant  as  a  function  of  temperature  allows  calculation  of  rates  for 
one  of  the  reactions  from  the  other.  At  lower  temperatures  the  heats  of  formation  of  all  four 
species  (two  reactants  and  two  products)  are  accurately  known,  thus  the  equilibrium  constant  can 
be  accurately  computed  at  these  temperatures.  Therefore,  the  good  agreement  between  the 


computed  rates  and  the  experimental  results  of  Demissy  and  Leselaux  for  reaction  (R5)  is 
inconsistent  with  the  ptx>rer  agreement  observed  between  theory  and  the  experimental  results  of 
Marshall  and  l  onti  jn  '-'  for  temperatures  below  660  K.  This  is  further  evidence  of  the 
inconsistence  of  the  lose  temperature  experimental  results  for  reaction  (R4). 

The  good  agreement  between  the  different  levels  of  theory  and  the  overall  good  agreement 
between  the  experimental  and  computed  rates  is  encouraging  for  appication  of  the  HF/BAC-MP4 
method  to  a  variety  of  reactions.  This  work  establishes  the  feasibility  of  this  approach  for 
providing  reliable  predictions  of  thermal  rate  constants  and  shows  the  utility  of  semiempirical 
corrections  to  ah  initio  electronic  structure  calculations  of  information  about  the  potential 
energy. surface.  In  particular,  the  HF/BAC-MP4  method  is  a  cost  effective  means  of  obtaining  this 
type  of  information.  Calculations  of  the  energy  and  its  first  and  second  derivatives  at  each  point 
along  the  minimum  energy  path  required  only  eight  minutes  of  epu  time  on  a  Cray  1  s  computer. 

To  obtain  rate  constants  numerically  converged  to  about  20 7c  at  300  K  requires  about  40  points 
along  the  reaction  path.  At  higher  temperatures  these  requirements  relax  so  that  at  l(XX)  K  only 
about  10  points  along  the  reaction  path  are  needed  for  the  same  type  of  convergence.  Thus,  for 
about  five  and  one-half  hours  of  Cray  Is  computer  time  rates  accurate  to  better  than  50T  can  be 
obtained  over  the  entire  temperature  range  from  200  to  24(X)  K.  The  calculations  at  the  MP2  and 
MP4  levels  become  much  more  costly  in  pan  because  some  or  all  of  the  analytically  derivatives  are 
not  available  and  must  be  computed  by  numerical  differentiation.  For  the  MP2  calculations  only 
analytic  first  derivatives  are  available  and  for  the  MP4  calculations  no  analytic  derivatives  are 
available.  Calculation  of  the  derivatives  by  finite  difference  greatly  increases  the  cost  of  the 
calculations.  The  MP2  and  MP4  calculations  require  over  three  and  seventeen  hours,  respectively, 
nf  Cray  Is  computer  time  per  point  along  the  minimum  energy  path.  The  reliabilitv  of  the  lower 
1  IF  level  of  theory  with  the  B  AC  correction  makes  this  efficient  approach  feasible. 


S.'s.'-’. 


CM* 


c.1.2  Application  to  HONO  —  UNO, 

As  a  further  test  of  the  BAC-MP4  method.  the  hydrogen  migration  reaction  from  trans- 
!  It  )\( )  to  UNO,  was  studied  at  several  levels  of  electronic  structure  theory.  Molecules  containing 

nitro  groups  present  a  challenge  to  ah  initio  methods  for  obtaining  accurate  representations  of  the 
potential  energv  surface.  The  utility  of  a  single  reference  wavefunction  has  been  questioned  for 
this  tvpe  of  s\  stem  and  it  has  been  suggest  that  a  multireference  wavefunction  such  as  used  in  the 
maitieoiitlguration  SCI-  t  MCSCF)  method  will  be  necessary  to  obtain  reasonable  estimates  of  the 
reaction  energetics.  Therefore,  this  type  of  reaction  represents  a  stringent  test  of  the  BAC-MP4 
method  vs  Inch  is  based  on  a  single  reference  I  IF  wavefunction. 

In  this  section  we  present  a  comparison  of  BAC-MP4  and  MCSCF  results  for  reaction 
'  R"7 1.  Two  different  MCSCF  calculations  are  considered,  both  using  double-zeta  plus  polarization 
!  D/P)  basis  sets.  A  DZP  basis  set  is  comparable  to  the  6-31G*  basis  set  used  in  the  HF 
optimizations  of  the  BAC-MP4  calculations  and  the  following  comparisons  are  a  test  primarily  of 
the  methixis  for  including  electron  correlation  and  not  a  test  of  basis  set  effects.  The  first  MCSCF 
calculation  treats  the  cr  and  rt  spaces  independently.  The  CT-space  correlation  includes  selected 
single  and  double  excitations  from  a  HF  reference  and  correlation  is  included  for  the  the  complete 
active  n  space;  this  is  denoted  MOa/tti.  The  second  MCSCF  calculation  considers  all  double 
exvitatiotis  from  the  seven  highest  occupied  orbitals  into  the  seven  lowest  unoccupied  orbitals;  this 
is  denoted  MCI 4. 

Tor  the  BAC-MP4  method  to  be  useful,  the  location  of  the  critical  geometries  (equilibrium 
e.i'metri'-s  and  saddle  points i  and  minimum  energy  path  for  the  reaction  must  be  accurately 
mpr  <d  cd  In  Injure  0  the  1  IF  and  MCio/tti  minimum  energy  paths  are  compared.  During  the 
.  ■c.re  reaction  all  tour  .it* >ms  he  in  a  plane;  figure  6  show s  the  motion  in  this  plane  of  the  1 1  and 
•.>.,)(  )  ,iii  mis  i  cl. it.  vo  ti  >  living  the  N  atom  at  the  origin  The  exeel  lent  agreement  of  the  two  paths 
ates  that  trie  B  \(  ’  MP4  methixl  has  passed  the  first  test;  the  HI  method  is  capable  of  placing 


F  igure  6.  Minimum  energy  path  for  the  unimolecular  reaction  MONO  — »  UNO,  computed  using  information  about  the 

potential  energy  surface  obtained  from  MCSCF  calculations  (solid  curves)  and  BAC-MP4  calculations  (•)  During  the 
entire  reaction  the  four  atoms  remain  in  a  plane.  The  motion  of  the  oxygen  and  hydrogen  atoms  in  this  plane,  relative  to 
the  nitrogen  atom  at  the  origin,  is  depicted. 


In  the  BAC-N1P4  methcxi  the  frequencies  are  also  calculated  at  the  HF  level  and  it  is 
interesting  to  see  how  well  they  compare  with  the  MCSCF  ones.  A  comparison  of  the  frequencies 
from  the  two  MCSCF  calculations  and  the  BAC-MP4  calculations,  as  well  as  comparison  to 
experiment40  is  presented  in  table  4. 


Table  4.  Frequencies  (in  cm'1)  for  trans  HONO 


Comparison  between  individual  modes  shows  some  fairly  large  discrepancies  but  the  crucial 
quantity  in  uetemiining  the  thermal  rate  is  the  total  zero  point  energy  given  at  the  bottom  of  each 
column.  Because  of  the  separation  of  the  a  and  k  spaces  in  the  MC(o/ti)  calculation  the  frequency 
of  the  lowest  mode,  which  is  an  out-of-plane  bend,  is  not  accurately  reproduced.  Otherwise  this 
MCSCF-  calculation  includes  the  most  correlation  and  gives  the  best  agreement  with  the 
experimental  zero  point  energy.  The  MC14  and  BAC-MP4  frequencies  agree  well  and  the  zero 
point  energy  for  these  methods  is  overestimated  by  about  189F  This  is  higher  than  the  usual  12 7c 
i  >\  eresiimate  expected  from  I  IF  frequencies  but  of  sufficient  accuracy  for  applying  the  BAC 
ci  irrcction.  The  theoretical  frequencies  at  the  saddle  point  are  presented  in  table  5.  In  this  case  the 
B  \0-.MP4  frequencies  agree  better  with  the  MC(a/7t)  ones  than  with  the  MC14  frequencies  but  the 


Table  5.  Frequencies  (in  cm'1)  for  the  migration  transition  state  HONO  — >  HN02. 


MC(o/7t) 

MCI  4 

B  AC- MP4 

2345i 

2385i 

2424i 

319 

436 

394 

341 

711 

709 

628 

1344 

1360 

765 

1680 

1686 

1179 

2629 

2434 

3234 

3400 

3292 

The  reaction  energetics  computed  from  the  BAC-MP4  and  MC(o/7t)  methods  are  presented 
in  figure  7  as  a  comparison  of  the  ground-state  adiabatic  potential  curves.  The  MCSCF  calculation 
has  errors  arising  from  basis  set  limitations  and  from  neglect  of  some  of  the  electron  correlation. 
The  uncertainties  in  the  energetics  obtained  from  this  method  are  estimated  to  be  about  10  kcal/mol. 
The  B  AC- MP4  gives  the  best  estimate  of  the  adiabatic  barrier  because  it  uses  a  bond  additivity 
correction  to  compensate  for  the  errors  caused  the  small  basis  set  and  the  limited  amount  of  electron 
correlation  included.  The  agreement  in  the  barrier  heights  may  be  fortuitous  (they  differ  by  less 
than  5  kcal/mol)  but  the  general  shape  of  the  two  curves  are  qualitatively  similar.  This  indicates 
that  the  BAC-MP4  method  is  not  giving  unphvsical  behavior  as  a  result  of  the  bond  additivity 
correction  over-  or  under-compensating  for  gross  errors  in  the  MP4  energies.  The  BAC-MP4 
method  gives  physically  reasonable  results  even  for  systems  containing  nitro  groups. 

Another  encouraging  feature  of  the  BAC-MP4  calculations  is  the  computational  effort 
required.  The  calculation  of  the  energy  and  first  and  second  derivative  at  each  point  along  the 
minimum  energy  required  less  than  20  minutes  of  Cray  Is  computer  time.  This  system  is  much 
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Figure  7.  Ground-state  adiabatic  potential  curves  as  a  function  of  reaction  coordinate  s  for  the 
unimoleeular  reaction  MONO  — >  HNO-i  computed  using  information  about  the  potential  energy 
surface  obtained  from  MCSCF  calculations  (dashed  curve)  and  BAC-MP4  calculations  (solid 
curve) 
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larger  than  the  H  +  NH-  reaction  considered  (24  electrons  as  compared  to  1 1)  but  the 


computational  effort  is  small  enough  to  allow  the  entire  reaction  path  from  -1.0  to  1.0  a  to  be 


mapped  out  in  less  than  eight  hours  of  Cray  Is  computer  time.  Thus,  we  are  encouraged  to  treat 


even  larger  polyatomic  systems  which  should  require  only  a  moderate  computational  effort. 


3.2  Systematic  studies  of  the  H  +  RNO-,  reactions  for  R  =  H,  CHV  OH,  NH-,,  and  CH-NH 


The  systematic  study  of  a  homologous  set  of  reactions  w  ith  different  substituents  allows 


for  an  assessment  of  the  effect  of  the  substituents  on  reaction  dynamics  and  is  the  beginning  of  any 


empirical  model  based  upon  group  additivities.  As  a  first  step,  such  a  study  has  been  performed 


on  reactions  (Rl)  and  (R2)  using  conventional  transition  state  theory  with  HF/BAC-MP4 


information  about  the  potential  energy  surface.  More  rigorous  studies  in  the  future  will  use 


variational  transition  state  theory.  For  reactions  (R3)  with  all  substituents  except  OH  the  HF/BAC- 


MP4  method  predicts  that  the  saddle  point  energy  is  below  the  asymptotic  energy  for  the  products. 


This  is  unphysical  and  requires  the  use  of  variational  transition  state  theory  for  a  consistent 


evaluation  of  the  thermal  rate  constants.  This  type  of  study  has  been  carried  out  for  one  of  the 


substituents  R=NH->  and  allows  a  complete  description  of  the  reaction  system  for  all  three  reactions 


(R1MR3)  for  the  substituents  H,  CHV  OH,  NIU,  and  CH,NH. 


3.2.1  H  +  R\0->  — >  RN( 


The  enthalpy  of  reaction  at  0  K,  A  H(),  is  the  change  in  the  ground-state  adiabatic  potential 


;urve  in  going  from  reactants  to  products.  This  gives  an  estimate  of  the  heat  release  (for  a  negative 


enthalpy)  upon  reaction  .  The  free  energy  of  reaction  for  a  given  temperature  T,  A  G-p  is  a  better 


measure  of  the  thermodynamic  driving  force  behind  the  reaction  since  it  includes  both  the  entropic 


and  enthalpic  terms.  Negatives  values  for  the  free  energy  of  reaction  indicate  that  the  products  are 


more  viable  than  reactants.  The  enthalpy  of  activation  at  0  K.  A  ll,  is  the  change  in  the  ground- 


o.ate  adiabatic  potential  curve  in  going  from  reactants  to  the  saddle  point  and  is  a  measure  of  the 


energetic  barrier  to  the  reaction.  Within  the  approximations  of  conventional  transition  state  theory 


I 
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the  rate  a;  which  the  reaction  pnveeds  for  a  given  temperature  is  governed  by  the  free  energy  of 
.tc;:\  ation.  A  G  j .  which  includes  the  energetic  factor  AtH()  and  the  entropic  terms.  The  enthalpies 

tree  energies  of  reaction,  the  enthalpies  free  energies  of  activation,  and  the  rate  constant  at  1  (XX)  K 
>r  reaction  <  R 1 )  is  summarized  in  table  6  for  all  five  substituents. 


Table  o. 

Reaction  energetics  and  rates  for  the  reaction  H  +  RNO-, 

P 

-4  RN( 

Y)h 

R 

V'n 

ArG-r 

AtGT 

k+(T=1000  K) 

1  kcal/mol  i 

(kcal/mol) 

(kcal/mol) 

(kcal/mol) 

(s'* 1) 

11 

-IS. 4 

X.7 

2X.3 

1.9x1  O'1- 

Cl  I; 

-41.1 

-17.4 

7.X 

31.1 

4.6x10' 13 

OH 

-35.0 

-16.2 

9.3 

30.0 

8 . 1  x  I  O' 1 3 

NIG 

32. ^ 

-12. s 

9.9 

29.X 

4.6x1  O'1’ 

Cl  G,  Nil 

->  -S  ■> 

-13.X 

9.0 

30.5 

6.3x1  O' 13 

it  is  clear  that  the  addition  of  a  hydrogen  atom  to  the  oxygen  end  of  RNO->  is  fairly  insensitive  to 

the  substituent  attached  to  the  nitrogen.  Although  the  enthalpy  of  reaction  changes  by  as  much  as  X 
Ks.tfmol  in  changing  the  substituent  from  Cl  I  ^  to  N1R,  the  enthalpy  of  reaction  changes  by  less 

than  2.2  heal  mol  for  all  five  substituents.  Similarly  changes  in  the  free  energies  of  reaction  are  as 
large  as  5,o  heal/mol  but  the  free  energy  of  activation  changes  by  less  than  3  kcal/mol.  The  rates 
change  at  most  by  about  a  factor  of  4  and  the  Arrhenius  parameters  are  ail  very  similar  at  1  (XX)  K  - 
the  actuation  energies  are  1 1.2.  10.2.  1  1.9,  1  1.3  and  1 1.3  kcal/mol  for  the  substituents  H.  CH^, 
OH.  NIG.  and  Cl  GN1 1,  respectively  and  the  base  ten  logarithms  of  the  A  factors  are  -9.3.  - 10. 1 ,  - 
0.5.  9  o.  ,md  -9.T  respectively. 

1  he  present  cals  illations  are  presented  as  a  qualitative  guide  to  substituent  effects.  A  more 
cons. stent  treatment  of  the  reaction  dynamics  would  involve  a  variational  transition  state  theory 
approach  including  tunneling  effects.  Since  these  reactions  involve  primarily  the  motion  of  a 
hsdrocen  atom  and  the  barriers  for  reaction  are  high  (ranging  from  S  to  10  kcal/mol),  tunneling  is 


expected  to  be  very  important  at  temperatures  below  about  6(X)  K.  At  1000  K,  where  the  current 


comparisons  are  being  made,  tunneling  should  enhance  the  rate  by  at  most  50%  and  the  qualitative 
comparison  made  above  should  be  valid.  For  meaningful  comparisons  at  temperature  below  600 
K  tunneling  should  be  treated. 

3.2.2  RN'(  -^R  +  HONO 
Ndh 

A  similar  comparison  of  reaction  energetics  and  rates  for  reaction  (R2)  is  provided  in  table 
7 . 
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Table  7.  Reaction  energetics  and  rates  for  the  reaction  RN(  — >  R  +  MONO 

XDH 


R 

V*o 

ArGT 

AtH(), 

AiGy 

k+(T=1000  K) 

(kcal/mol) 

(kcal/mol) 

(kcal/mol) 

(kcal/mol) 

(s'1) 

11 

28.9 

5.7 

34.4 

34.3 

6.6x1 05 

CH, 

20.2 

-1 1.6 

28.3 

25.7 

5.0x1 07 

OH 

3.0 

-25.4 

6.0 

6.2 

9.2xlOn 

NHi 

4.2 

-28.1 

10.7 

10.4 

1 . 1  x  1 0 1 1 

CH;NH 

4.6 

-31.3 

10.9 

11.7 

5.7x1 01() 

The  rates  reported  for  the  unimolecular  reaction  are  for  the  high  pressure  limit.  These  reactions 
inv  olve  the  breaking  of  the  R-N  bond  and  the  thermal  rates  are  very'  sensitive  to  the  substituent  R. 
There  is  a  weak  correlation  between  the  free  energy  of  reaction  and  the  free  energy  of  activation 
anil  thereby  the  reaction  rate.  More  important  though  is  the  comparison  between  the  energetics  and 
rates  for  the  NI I->  and  Cl  RNl 1  substituents.  In  these  cases  a  NO  bond  is  being  broken  in  the 

course  of  the  reaction  and  the  energetics  and  rates  tire  very  similar. 


/O  P 

3.:..'.  N1KN(  ->  Ml,  +  IlONO  and  NH->N(  ->NH,NO  +  OH 

'  V)H  -  XDH 

An  important  aspect  of  the  reaction  system  (R1HR3)  is  the  branching  ratio  between  the 
two  uniniolecular  reactions  (R2)  and  (R3F  This  has  been  studied  in  more  detail  for  the  substituent 

P 

M  1,.  Rate  constants  tor  the  reaction  NH,N(  — >  \H,NO  +  OH  were  computed  using 

variational  transition  state  theory  with  semiclassical  adiabatic  ground-state  transmission  coefficients 
using  I IF/BAC-MP4  information  about  the  potential  energy  surface.  The  potential  along  the 
minimum  energy  path  and  the  ground-state  adiabatic  potential  curve  are  shown  in  figure  8.  Note 
that  the  value  of  the  potential  at  the  HF  saddle  point  location  (  indicate  by  s=0)  is  much  lower  than 
its  maximum  value  which  occurs  for  s  at  about  0.25  a  .  This  same  type  of  behavior  is  seen  in  the 

ground-state  adiabatic  potential  curve  and  in  the  free  energy  of  activation  curves  show  in  figure  9 
for  300  and  1000  K.  As  a  result,  the  conventional  transition  theory  state  rates  are  much  higher 
than  the  ones  obtained  from  the  improved  canonical  variational  theory  as  seen  in  table  8. 


Table  8. 

Rate  constants  (in  units  of  s'1 

)  for  the  two  reactions 

NH,> 

/° 

«'(  -4  NH 

XDI I 

2NO  +  OH  NH,n 

,0 

.(  — »  NH,  +  MONO 

ull 

T.K 

k'(T) 

kIcvr(T) 

|.ICVT/SCSAG(T) 

kf(T> 

300 

9.4x10° 

3.5x1 O6 

5.3x1 O6 

l.lxl  0^ 

400 

6.9x10'° 

2.4x10s 

3.0x10s 

1.4xl()7 

600 

5.6x10" 

1.9x10'° 

2.0x10'° 

1.9x10° 

1000 

3.4\10'“ 

7.4x10' 1 

7.6x10" 

1.1x10" 

1500 

S3K1012 

5. Ox 1012 

5.  lx  10 12 

9.4x10' ' 

2400 

1.9x10" 

1.7x10' ' 

1.7x10" 

5.2x1  ()' 2 

:£v:v:v:>:v:v:^ 
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Figure  s.  Potential  energy  along  the  minimum  energy  path  [  Vyj|.-p(s)|  and  ground-state  adiabatic 

potential  [  VAi  si|  as  a  function  of  reaction  coordinate  s  computed  using  the  I IF/BAC-N1P4  method 
for  the  reaction  MIA'iOiOH  -4  NH-.NO  +  on. 
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Tunneling  is  relatively  unimportant  for  this  reaction  and  the  rate  obtained  by  the  1CVT/SCSAG 


method  are  onlv  sliuhtlv  higher  than  those  with  the  tunneling  correction  factor. 

'  ‘  ^  /O 

I  he  rates  tor  the  reaction  NTHNX  — >  NI I-,  +  HONO  are  also  presented  in  table  8  for 

comparison.  The  value  of  the  ground-state  adiabatic  barrier  above  the  reactant  value  is  about  10.7 
keal/mol  at  the  saddle  point  for  the  HONO  channel  and  about  9.1  keal/mol  at  its  maximum  for  the 
01 1  channel.  A  consistent  comparison  would  treat  both  branches  of  the  reaction  using  variational 
transition  state  theory;  however,  if  this  is  done,  the  value  of  the  adiabatic  barrier  maximum  will 
increase  for  the  HONO  channel  making  it  even  more  dynamically  unfavorable.  Since  the  rates  for 
the  OH  channel  tire  already  higher,  increasing  the  barrier  for  the  HONO  channel  will  only  increase 
the  difference.  Therefore,  we  can  predict  that  the  OH  channel  is  favored  over  that  for  HONO 
formation  at  all  temperatures  from  300  to  2400  K. 

3.3  Thermochemical  kinetic  analysis 

Appendix  B  gives  a  more  indepth  discussion  of  the  results  of  the  thermochemical  analysis. 
In  this  section  we  review  those  results  and  present  additional  analysis  for  reaction  (R4)  and  (R6). 
The  new  thermochemical  kinetic  analysis  presented  in  detail  in  appendix  B  includes  important 
dynamical  effects  that  will  allow  reliable  estimates  of  Arrhenius  parameters.  As  a  test  of  the 
importance  of  these  dynamical  effects  and  as  a  guide  to  estimating  the  nonsubstantial  contributions 
the  thermochemical  kinetic  analysis  is  performed  on  seven  reactions.  Five  of  the  reactions  have 
been  treated  earlier  by  variational  transition  state  theory10,11,13'15  and  provide  rate  constants  in 
g<x>d  agreement  with  experiment,14,41  and  the  other  two  reactions  are  reaction  (R4)  and  (R6) 

already  discussed  in  this  report.  The  theoretical  rate  constants  are  analyzed  to  obtain  values  for 
A n S  AnIIl|.  A,.S't  and  AJl'-j..  The  five  new  reactions  considered  are 


Q(  P)  +  — >  OH  +  H 


OH  +  H,  ->  11,0  +  H 


•or  those  five  reactions  global  ab  initio  potential  energy  surfaces  were  used  in  contrast  to  the 


reaction  path  I  lamiltonian  information  used  for  reactions  (R4)  and  (R6). 


The  activation  parameters  for  the  seven  reactions  are  present  in  table  9  and  are  given  in 


terms  of  the  entropy  and  enthalpy  of  activation  at  300  K,  and  the  heat  capacity  of  activation  from 


200  to  1500  K.  Some  general  trends  are  seen  in  table  9.  First,  using  the  terminologv  discussed  in 


section  2.3,  the  substantial  contributions  are  large  than  the  nonsubstantial  ones;  thus,  the  analysis 


in  terms  of  the  structure  of  the  transition  state  gives  the  majority  of  the  contributions  to  the  total 


mthalpv  and  entropy  of  activation.  However,  the  nonsubstantial  contributions  to  both  the  enthalpy 


and  entropy  of  activation  are  quite  large  and  need  to  be  included  for  accurate  predictions  of  the  rate 


constants  The  values  of  the  substantial  enthalpy  of  activation  at  300  K  correlate  with  the  barrier 
heights  for  the  reactions,  that  is.  OH  +  H-,  has  the  smallest  barrier  and  H  +  NH^  has  the  largest. 


For  the  nondeuterated  reactions  the  nonsubstantial  contributions  tend  to  be  larger  for  the  svstems 


v.  ith  the  larger  barriers,  i.e.,  tunneling  is  more  important  for  the  reactions  with  the  higher  barriers. 


Compared  to  the  hydrogen  cases,  the  deuterated  reactions  have  smaller  nonsubstantial 


contributions  due  to  the  lesser  importance  of  tunneling. 


The  substantial  contributions  to  the  heat  capacities  of  activation  are  negative,  in  general. 


which  implies  that  the  substantial  enthalpies  and  entropies  decrease  with  increasing  temperature.  In 


the  present  example,  the  nonsubstantial  contributions  to  the  heat  capacity  of  activation  go  to  zero  at 


high  temperatures.  They  are  all  positive  and  at  low'  temperatures  are  much  larger  than  the 


mbstantut!  ones  causing  large  changes  in  the  enthalpies  and  entropies  at  low  temperatures. 


The  substantial  contributions  to  the  activation  parameters  are  obtained  from  information 


alxmt  the  potential  at  a  single  point  along  the  minimum  energy  path.  This  enables  the  substantial 


c.tnbutions  to  lx-  estimated  from  existing  bond  additivity  or  group  additivity  relationships/'  The 
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uonsubstantial  contributions  arc  obtained  from  more  extended  information  about  the  potential 
nacres  surface  and.  therefore,  are  harder  to  estimate.  All  the  reactions  studied  here  are  hydrogen 
•  or  deuteriums  atom  transfer  reactions  with  fairly  high  barriers.  l  or  these  reactions  the 
uonsubstantial  contributions  are  predominantly  from  quantum  mechanical  tunneling  and  the  effect 
of  \  ariationally  optimizing  the  location  of  the  transition  state  dividing  surface  is  small.  For 
reactions  with  smaller  barriers  we  expect  the  effect  of  quantum  mechanical  tunneling  to  be  smaller, 
mu  for  those  reactions  the  effects  of  varuttionally  locating  the  dividing  surface  w  ill  become  more 
important  especially  at  higher  temperatures.4-  It  is  therefore  expected  that  the  uonsubstantial 
c.  'mr;  but  ions  will  be  significant  for  a  w  ide  range  of  gas-phase  chemical  reactions  and  new'  methods 
are  needed  tor  estimating  these  quantities. 

4.  Conclusions 

The  main  objective  of  this  six-month  Phase  I  research  program  was  to  establish  the 
feasibility  of  using  theoretical  methods  for  describing  the  kinetics  of  chemical  reactions  important 
in  nitrocubane  decomposition.  The  accuracy  of  the  theoretical  methods  has  been  established  by 
comparison  of  different  levels  of  theory  for  two  chemical  reactions,  H  +  NH-  — >  H->  +  NH-.  and 

I  iONO  — >  UNO-)  and  in  comparison  of  the  theoretical  rate  constants  w  ith  experimental  ones  for  the 

I I  i-  M  I;  — » II-,  +  NIC  reaction.  These  calculations  indicate  that  the  underlying  ab  initio 
electronic  structure  calculations  provide  a  good  zeroth- order  description  of  the  potential  energy 
surfaces.  The  empirical  bond  additivity  correction  compensates  for  known  errors  in  these 
potentials  and  the  comparisons  w  ith  experiment  indicate  that  these  correction  are  of  quantitative 
accuracy. 

For  these  theoretical  methods  to  be  of  practical  use  they  must  be  capable  of  treating  systems 
w  .-.re  larger  than  those  presented  here.  The  only  inherent  limitation  in  these  methods  is  the 
,  ’ ■ :  ;  .t.itiotial  effort  of  the  ah  initio  electronic  structure  calculations  which  increases  with 
.mm,:-.:.:  s;/e  of  the  reaction  svstem.  A  requirement  for  the  utility  of  this  method  is  that  the  size 
me  ah  :>m„>  electronic  structure  calculations  is  kept  manageable.  Although  this  leads  to  known 
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errors  in  the  potential  energy  surfaces,  the  bond  additivity  correction  provides  a  computationally 
efficient  method  of  compensating  for  these  errors.  The  relatively  small  reaction  systems  studied 
here  were  treated  with  a  very  modest  expenditure  of  computational  power,  e.g.,  the  entire 
calculation  of  the  rate  constants  from  2(X)  to  2400  K.  for  the  H  +  NH^  — >  H->  +  NHt  reaction  were 

carried  out  in  several  hours  of  Cray  Is  computer  time.  The  prospects  of  treating  much  larger 
systems  are  excellent. 

Looking  beyond  the  Phase  I  objectives  towards  the  goal  of  providing  a  comprehensive 
treatment  of  propellant  decomposition,  methods  are  needed  for  treating  the  kinetics  of  very’  large 
poly  atomic  systems.  Accurate  calculations  of  rates  for  reactions  involving  very'  large  polyatomic 
molecules  is  beyond  the  current  computational  capabilities  of  the  electronic  structure  methods  even 
w  ith  the  use  of  empirical  corrections.  However,  the  kinetics  of  these  systems  can  still  be 
addressed  using  approximate  methods  for  estimating  kinetic  parameters  based  on  bond  additivity 
and  group  additivity  relationships.  This  type  of  approach  requires  the  development  of  a  database 
of  information  about  substituent  effects  on  kinetic  parameters.  As  a  first  step  in  this  direction,  a 
s\  stematic  study  was  carried  out  for  hydrogen  attack  on  simple  nitro  containing  compounds. 

Three  different  reaction  paths  were  considered  and  the  results  show  that  the  rates  for  one  of  them 
are  relatively  insensitive  to  substituent  effects,  whereas  the  rates  for  the  other  two  change  by 
several  orders  of  magnitude  upon  interchange  of  substituents. 

Qualitative  information  about  substituent  effects  is  in  itself  valuable;  however,  to  be  of 
practical  use  for  providing  reliable  estimates  of  kinetic  parameters  for  large  polyatomic  systems, 
methods  of  quantifying  this  type  of  information  are  needed.  Benson's^*  thermochemical  kinetic 
analysts  of  chemical  reactions  provides  such  a  quantification  scheme.  Rates  are  based  upon 
activation  parameters  (the  enthalpy  and  entropy  of  activation)  which  are  obtained  from  structural 
properties  of  the  saddle  point.  These  can  often  be  reliably  estimated  using  bond  or  group  additivity 
relationships.  The  main  deficiency  in  this  approach  is  that  it  neglects  important  variational  and 
quantum  tunneling  effects.  A  new  approach  to  thermochemical  kinetic  analysis  has  been  presented 
in  tins  report  which  shows  the  great  importance  of  these  effects  and  provides  a  prescription  for 
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including  them  in  the  analysis.  The  results  reported  here  are  the  beginning  of  a  new  database  of 
kinetic  parameters  which  includes  all  the  important  dynamical  effects  needed  for  accurate 
calculations  of  the  rates.  Calculations  of  activation  parameters  for  smaller  systems  using  the 
semiempirical  methods  validated  in  this  report  will  be  used  to  build  up  this  database  to  allow 
treatment  of  svstems  involving  very  large  polyatomic  molecules. 
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Appendices  A. 

Theoretical  calculations  of  the  thermal  rate  constants  for  the  gas-phase 
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Theoretical  calculations  of  the  thermal  rate  constants  for  the  gas-phase 
chemical  reactions  11  +  NH^  H->  +  NH-,  and  D  +  ND^  ^  D->  +  ND-, 
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Abstract.  Rate  constants  for  the  title  reactions  are  computed  using  variational  transition  state  theory 
with  semiclassical  ground-state  adiabatic  transmission  coefficients  for  the  temperature  range  from 
200  to  2400  K.  The  rates  are  computed  from  limited  information  about  the  potential  energy  surface 
along  the  minimum  energy  path  as  parameters  of  the  reaction  path  Hamiltonian.  The  potential 
information  is  obtained  from  ab  initio  electronic  structure  calculations  with  an  empirical  bond 
additivity  correction.  The  accuracy  of  this  semiempirical  technique  for  obtaining  the  potential 
information  is  tested  by  using  increasingly  higher  quality  ab  initio  electronic  structure  calculations 
before  applying  the  Kind  additivity  correction.  For  the  reactions  H  +  NHj  H-,  +  NH2  and  D  + 
ND,  — >  D-,  +  ND-,  the  ultimate  test  is  given  by  comparison  with  recent  experimental  results.  The 
agreement  is  very  good  except  for  the  very  low  temperature  results  for  H  +  NH^  — >  H-,  +  NH2; 
the  excellent  agreement  for  the  reverse  reaction  in  the  same  temperature  range  indicates  an 
inconsistency  in  these  experimental  results. 


. 


!.  Introduction. 

I  lydrogen  atoms  abstraction  reactions  are  ubiquitous  in  the  combustion  of  hydrogen 
containing  compounds.  As  an  example,  the  gas-phase  chemical  reaction 

II  ~  MI; II,  +  Ml,  (Rl) 

is  important  in  the  tiiennal  decomposition,  pyrolysis,  and  combustion  of  ammonia.  The  rate  of 
reaction  { R 1  >  has  recently  been  measured  using  three  different  experimental  techniques;1-’  these 
experiments  combine  to  yield  rate  data  from  500  K.  -  where  the  slowness  of  the  reaction  tests  the 
limits  of  the  sensitivity  of  the  experimental  apparatus  -  up  to  1777  K.  These  experiments  are  in 
good  agreement  w  ith  each  other  in  temperature  regions  where  they  overlap,  but  overall  they 
disagree  with  older  experiments4-'1  that  indicated  the  reaction  was  much  slower. 

Compared  to  reaction  (Rl ).  its  reverse 

11,  -  NIC  -»  H  +  NH;  (R2) 

has  been  the  subject  of  fewer  direct  measurements'-6’'  although  it  is  also  important  in  ammonia 
combustion.  The  most  recent  measurement  of  the  rate  of  reaction  (R2)  extends  over  the 
temperature  range  from  673  to  1003  K.  using  a  single  experimental  technique4  and  combined  with 
the  earlier  study6  provide  reactions  rates  from  400  K  up  to  1003  K.  By  also  studying  (Rl)  in  the 
same  experimental  apparatus,  new  estimates  of  the  thermochemistry  of  the  Nil-,  radical  have  also 

been  obtained.  This  type  of  information  is  useful  in  evaluating  the  equilibrium  constant  and  in 
relating  the  rates  of  reactions  (Rl )  and  (R2)  by  detailed  balance. 

The  deuterated  analog  of  reaction  ( R 1 ) 

D  +  ND;  — »  D,  +  NT>2  (R3) 

lias  also  been  studied  using  a  single  experimental  technique  over  the  temperature  range  from  590  to 
1  220  KT  The  experimental  rates  of  reactions  (Rl)  and  (R3)  have  been  compared  with  theoretical 
ones  bawd  upon  conventional  transition  state  theory.^  These  calculations  utilize  information  about 
the  potential  energy  surface  obtained  from  electronic  structure  calculations1*-8  using  a  bond 
additivity  correction  t  the  BAC-Mi'4  method).101"'  The  agreement  of  the  theory  and  experiment 
indicates  that  the  calculated  barrier  height  of  about  IS  keal/mol  (or  about  16  keal/mol  when 


corrected  tor  changes  in  zero  point  energy)  is  reasonable.2  Curvature  in  the  Arrhenius  plot  of  the 
rate  for  reaction  (R 1 )  at  low  temperatures  has  been  attributed  to  quantum  mechanical  tunneling,  and 
attempts  have  been  made  to  substantiate  this  claim  using  simple  tunneling  models. 

The  recent  experimental  work  on  the  reactions  (R 1  )-(R3)  provides  an  excellent  testing 
ground  for  theoretical  methods  for  predicting  gas-phase  reaction  rates  since  all  three  reactions 
involve  the  same  potential  energy  surface.  One  of  the  purposes  of  the  present  paper  is  to  provide  a 
more  rigorous  treatment  of  the  dynamics  controlling  rates  of  reactions  (Rl)  and  (R3)  and  to  also 
provide  the  first  theoretical  estimates  of  the  rate  for  reaction  (R2)  and  predictions  for  its  deuterated 
anttlog 

D-i  +  ND-)  — ?  D  +  ND^  (R4) 

Previous  calculations  on  the  reaction  (Rl)  and  (R3)  employed  conventional  transition  state  theory 
in  which  the  dynamical  bottleneck  is  assumed  to  occur  at  the  saddle  point  for  the  reaction.  As  will 
be  show  n  later,  because  of  the  nature  of  the  electronic  structure  information  about  the  potential 
energy  surface  utilized  in  these  calculations,  the  use  of  conventional  transition  state  theory  leads  to 
inconsistencies  in  the  treatment;  but  by  basing  the  calculations  upon  variational  transition  state 
theory1-'22  these  inconsistencies  are  no  longer  a  problem.  Another  advantage  of  using  variational 
transition  state  theory  is  that  it  allows  for  a  consistent  treatment  of  quantum  tunneling  effects 
which  are  important  at  low  temperatures.  A  second  goal  of  this  work  is  to  reexamine  the 
importance  of  quantum  mechanical  tunneling  in  these  reactions  utilizing  consistent  semiclassical 
adiabatic  transmission  coefficients  and  extending  the  rate  data  to  lower  temperatures. 

Compared  to  other  dynamical  methods  for  obtaining  reaction  rates,  variational  transition 
state  theory  has  the  practical  advantage  of  not  requiring  a  global  representation  of  the  potential 
energy  surface.  From  knowledge  of  the  potential  in  a  region  around  the  minimum  energy  path 
connecting  reactants  to  products,  accurate  estimates  of  the  reaction  rate  can  be  made.  This  type  of 
potential  information  can  be  directly  obtained  from  modern  electronic  structure  methods  as  the 


energy  and  derivatives  of  the  energy.  Previous  conventional  transition  state  theory  calculations  of 
the  rates  for  reactions  (RH  and  (R2)  have  utilized  potential  information  obtained  by  the  BAC-MP4 
method.10' 14  One  of  the  most  important  purposes  of  this  paper  is  to  critically  evaluate  the  method 
used  to  obtain  this  potential  energy  surface  information. 

2.  Theory 


2.1.  Variational  transition  state  theory 

Variational  transition  state  theory  with  semiclassical  transmission  coefficients  have  been 
extensively  review  in  the  literature.1-'  --  In  the  present  work  we  briefly  review  the  method  with  an 
emphasis  on  the  potential  energy  surface  information  needed  for  the  calculations. 

Conventional  transition  state  theory4  reduces  the  calculation  of  the  rate  constant  to  one  of 
quasiequilibrium  statistical  mechanics;  the  equilibrium  rate  is  approximated  as  the  equilibrium  flux 
headed  towards  reactants  through  a  dividing  surface  located  at  the  saddle  point.  With  this 
approximation,  the  rate  constant  kdT)  for  temperature  T  takes  on  the  simple  familiar  form 


k  ( T )  = 


kBT 
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Q  ( T ) 
.  K 
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where  ku  is  Boltzmann's  constant,  h  is  Planck's  constant.  Q'(T)  is  the  partition  function  for  the 

t  )  ,  4-  , 

bound  degrees  of  freedom  at  the  saddle  point,  Q>  (T)  is  the  reactants  partition  function,  and  V'  is 
the  value  of  the  potential  at  the  saddle  point.  Thus  conventional  transition  state  theory  requires 
information  about  the  potential  energy  surface  only  in  the  saddle  point  and  reactant  regions.  If  the 
partition  functions  are  computed  using  a  harmonic  approximation  then  the  matrix  of  second 
derivatives  (the  Hessian  matrix)  suffices. 

In  v  ariational  transition  state  theory  the  dividing  surface  is  viewed  as  a  tentative  dynamical 
bottleneck  to  flux,  and  the  best  bottleneck  (the  dividing  surface  allowing  the  least  flow  of  flux)  is 
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located  by  a  \anational  procedure.  The  generalized  expression  for  the  thermal  rate  constant  is 
gixen  as  a  function  of  the  location  s  of  the  dynamical  bottleneck  along  the  reaction  coordinate 
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where  QCl  1  tT.s)  is  the  generalized  partition  function  for  the  bound  degrees  of  freedom  orthogonal 
to  the  reaction  path  at  s.  and  V^pls)  is  the  value  of  the  potential  along  the  reaction  path  at  s.  One 

version  of  variational  transition  state  theory,  the  canonical  variational  theory  (CVT),-4  results  from 
minimizing  eq.  (2)  with  respect  to  s 

,  CVT  .  GT,  _ 

k  (  f)  =  nnn  k  (T,  s)  ^ 

The  improved  canonical  variational  theory  tICVT)-4  also  variationally  optimizes  the  location  of  the 
transition  state  dividing  surface  for  a  given  temperature,  but  provides  an  improved  treatment  of 
threshold  energies  by  using  an  ensemble  which  removes  energies  below  the  ground-state  adiabatic 
threshold.  To  compute  the  rate  constant  using  either  the  canonical  or  improved  canonical 
variational  theory  more  information  about  the  potential  energy  surface  is  required  than  for  a 
conventional  transition  state  theory  calculation;  information  about  the  potential  in  a  region  around 
the  reaction  path  is  required. 

The  reaction  path  is  defined  as  the  minimum  energy  path  connecting  the  saddle  point  with 
birth  the  reactant  and  product  regions.  The  minimum  energy  path  is  located  by  following  the  path 
of  steeper!  decents  in  both  directions  from  the  saddle  point  in  a  mass  weighted  coordinate  system 
such  that  each  degree  of  freedom  has  the  same  effective  mass  in  the  kinetic  energy  expression.  The 
first  step  from  the  saddle  point  is  taken  along  the  eigenvector  of  the  negative  eigenvalue  of  the 
1  lessian.  and  subsequent  steps  are  taken  in  the  direction  of  the  negative  of  the  gradient  (the  first 
derivative  of  the  energy).  lor  a  harmonic  treatment  of  the  partition  functions,  the  matrix  of  second 
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derivatives  along  the  minimum  energy  path  will  sulfiee.  In  this  ease  the  potential  information 
needed  is  the  energy  and  its  first  and  second  derivatives  along  the  minimum  energv  path. 

In  the  above  expressions  for  the  rate  constant  the  hound  degrees  of  freedom  are  treated  quantum 
mechanical's  i:.e..  the  partition  functions  are  evaluated  quantum  mechanically)  hut  the  motion 
.hong  trie  reaction  coordinate  (the  degree  of  freedom  corresponding  to  the  negative  eigenvalue  of 
the  1  lessian  at  the  saddle  point)  is  treated  classically.  Quantum  mechanical  effects  on  the  reaction 
Cisudmate  motion  le.g..  quantum  mechanical  tunneling)  is  included  by  a  multiplicative  factor  -  the 
transmission  coefficient.  The  Wigner  correction  factor^  has  been  very  popular  for  estimating 
quantum  mechanical  tunneling  effects  because  it  requires  know  ledge  of  the  shape  of  the  classical 

iwir  t!  ic  saddle  point;  this  type  of  information  can  be  obtained  from  the  negative  eigenvalue 
of  the  I  Icssi.m  at  the  saddle  point.  I  nfortunately,  the  approximations  of  this  method  are  rarely 
v.b.d  for  chemical  reactions,  i.e..  when  tunneling  is  important  it  occurs  over  distances  much  longer 
than  that  represented  by  a  quadratic  representation  of  the  potential  near  the  saddle  point,  and  a  more 
global  description  of  tunneling  is  required. 

A  consistent  route  to  including  quantum  mechanical  effects  on  reaction  coordinate  motion  is 
pr- >\;dcd  in  the  v ibrationally  and  rotational!)  adiabatic  theory  of  reactions. When  reaction 
coordinate  m.  non  is  treated  classically  one  form  of  variational  transition  state  theory'  -  the 
microcanoiiical  variational  theory  yields  an  expression  for  the  thermal  rate  constant  which  is 
equivalent  to  that  obtained  from  the  adiabatic  theory  of  reactions  Thus  it  is  consistent  within 
v  .;riu&>::.d  transition  state  theory  to  treat  tunneling  as  occurring  through  the  one-mathematical- 
d. mats., \  .b:\itionally-n nationally  adiabatic  potential 


V  s.  (/.  i  =  V  , 
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use  index  of  the  quantum  numbers  for  the  bound  modes  and  t'V,  (s)  is  the 


level  for  si.uc  u  at  the  cenerah/ed  transition  state  located  at  s  alone  the  reaction  path. 


For  thermal  rate  constants  the  tunneling  is  approximated  using  only  the  ground  state  adiabatic 
potential  curve  (U=0)  which  is  denoted  V‘Vs)  The  adiabatic  approximation  is  made  in  a 
curvilinear  coordinate  system,  and,  although  the  potential  term  is  very  simple,  the  kinetic  energy 
term  is  complicated  by  factors  dependent  upon  the  curvature  of  the  reaction  path.  For  systems  in 
which  the  curvature  of  the  reaction  path  is  not  too  severe,  the  small-curvature  semiclassical 
adiabatic  ground  state  (SCSAG)  method  '’0  includes  the  effect  of  the  reaction-path  curvature  to 
induce  the  tunneling  path  to  'cut  the  corner'  and  shorten  the  tunneling  length.  Combining  the 
SCSAG  transmission  coefficient  with  the  improved  canonical  variational  theory  rate  constant  yields 

k!CVT/SCSAC(T)  =  icSCSAG(T)klCVT(T)  (4) 

To  construct  the  adiabatic  potential,  the  type  of  potential  information  required  is  identical  to  that 
needed  for  the  variational  transition  state  theory  calculation.  For  the  SCSAG  calculation  it  is  also 
necessary  to  know  the  curvature  of  the  reaction  path  which  can  be  obtained  from  the  first  and 
second  derivatives  of  the  potential  along  the  reaction  path.  Thus,  to  provide  a  consistent  and 
accurate  estimate  of  the  tunneling,  no  new  information  about  the  potential  energy  surface  is 
required. 

The  information  about  the  potential  energy  surface  is  used  in  the  POLYRATE  program  in 
the  form  of  parameters  of  the  reaction  path  Hamiltonian31  at  a  set  of  points  along  the  reaction  path. 
Because  of  the  expense  of  the  electronic  structure  calculations  the  information  about  the  potential  is 
often  limited,  and  the  parameters  of  the  reaction  path  Hamiltonian  are  interpolated  from  a  sparse 
grid  and  can  be  extrapolated  out  to  the  asymptotic  reactant  and  product  regions.15,32 

2.1.  Quantum  chemistry 

Despite  tiie  developments  in  quantum  chemistry'  in  the  past  decade,  it  remains  a  challenge  to 
calculated  accurate  potential  energy  surfaces  useful  for  determining  rates  of  activated  chemical 
reactions.  This  is  largely  due  to  the  fact  that  small  changes  in  the  barrier  height  translate  into  large 

(i 


changes  in  the  reaction  rate  (e.g.,  a  change  of  only  0.5  keal/mol  at  room  temperature  leads  to  a 
change  in  the  rate  of  over  a  factor  of  two).  Currently,  ab  initio  calculations  of  energetics  to  within 
chemical  accuracy  is  affordable  only  for  relatively  small  systems.  Even  so,  progress  has  been 
made  in  developing  methods  for  determining  reaction  energetics,  but  these  have  relied  on  empirical 
corrections  to  ab  initio  calculations. 

The  concept  of  a  molecule  as  a  collection  of  individual  bonds  between  atoms  is  central  to 
the  chemist's  view  of  molecular  structure  and  its  rearrangement  during  reactions.  The  order  of  a 
bond  (e.g.,  single  bonds  between  monovalent  atoms,  double  bonds  between  divalent  atoms,  etc.) 
has  been  used  in  establishing  empirical  rules  relating  bond  order  to  bond  lengths33  and  bond 
strengths.-  The  empirical  knowledge  that  the  dissociation  energy  for  a  given  bond  with  a  given 
bond  order  is,  in  general,  roughly  the  same  in  different  molecules  has  lead  to  the  chemist’s 
intuition  that  bond  strengths  tire  approximately  additive.  This  concept  of  bond  additivity  has  been 
successfully  applied  to  obtaining  estimates  of  the  thermochemistry  of  stable  species.35  Within  the 
last  several  years,  this  concept  has  been  used  in  a  more  quantitative  approach  (the  BAC-MP4 
method)  in  which  corrections  are  made  to  ab  initio  calculations.10'14  The  BAC-MP4  method  has 
been  successfully  applied  to  predicting  the  thermochemistry  of  a  wide  variety  of  stable  species  and 
transient  radical  species  that  occur  in  combustion  environments.  More  recently,  the  method  has 
been  extended  to  the  calculation  of  the  energetics  of  transition  states  for  chemical  reactions.1’^-36 

Inherent  in  the  BAC-MP4  method  is  the  assumption  that  geometries  and  frequencies  can  be 
obtained  to  a  good  approximation  from  a  Hartree-Fock  calculation  using  a  moderate  sized  basis. 
Therefore,  a  BAC-MP4  calculation  begins  with  a  Hartree-Fock  (HF)  optimization  to  the 
equilibrium  or  transition  state  geometry  (i.e.,  locating  the  geometry  at  which  the  gradient  vanishes) 
with  a  contracted  atomic  orbital  basis  which  for  the  present  application  is  denoted  6-31G*.37  The 
I  iessi.m  matrix  ,s  also  computed  at  the  I  IF  level  using  the  same  basis.  The  absolute  energy  fora 
given  geometry  is  much  harder  to  accurately  compute  and  requires  including  a  large  percentage  of 
the  electron  correlation  energy.  As  a  first  step  in  this  direction,  at  the  optimized  HF  geometry,  the 
energy  is  computed  by  full  fourth-order  Moller-PIesset  perturbation  theory3,3’39  (MP4-SDTQ), 
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including  the  effects  of  all  single,  double,  triple,  and  quadruple  excitations  in  the  contracted  atomic- 
orbital  basis  denoted  6-31G**.-7  Although  the  MP4  method  does  include  a  large  fraction  of  the 
electron  correlation  energy,  there  are  still  errors  due  to  basis  set  and  higher  order  correlation 
effects,  so  the  bond  additivity  correction  (BAC)  is  applied. 

The  BAC  correction  is  given  as  the  sum  of  contributions  over  all  bonds  in  the  molecule 


,  ,r  =  -  ( A  -  EMP4)  -  Z  Z  A.,  exp  (-a..  R  A  FI  f.. 
5AC  i  v  1  01  1  '  i  j  <  i  *J  v  u  ij)  kiti.  ij 
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where  AH^  is  experimental  atomic  heat  of  formation  of  atom  i,  E^4  is  the  atomic  energy 
computed  at  the  MP4  level  for  atom  i,  A-  and  a-  are  parameters  determining  the  strength  and 
range  of  the  correction,  R-  is  the  length  of  the  bond  between  atoms  i  and  j,  and  the  terms  fijk  damp 
out  the  bond  additivity  correction  when  atom  k  is  bonded  to  either  atom  i  or  j.  This  damping  factor 
takes  the  form 


1  ,jk  =  { 1  -  F  jF  k  «P  [  -  lk(  R  ik  -  R 0>]  }  { 1  -  F . F  k  txP  L“  2 “  jk(  R jk  ^  R  o)  J} 


where  the  parameters  Ft  are  dependent  only  on  the  type  of  atom  i.  A  bond  is  assumed  to  be  formed 
between  atoms  i  and  j  if  the  bond  length  is  less  than  some  cutoff  value 


R  .  <  1.  3  Ral  +  Ral 


where  the  parameter  R^  is  a  measure  of  the  atomic  radius  of  atom  i. 

The  parameters  of  the  BAC-MP4  method  are  fitted  to  a  database  of  experimental  heats  of 
formation  of  atomic  and  molecular  species  in  a  least  squares  sense.  First  the  parameters  AE,  tire 
fitted  to  reproduce  atomic  heats  of  formation,  then  the  parameters  A(|,  cq.,  and  F(  are  fitted  to 


reproduce  diatomic  and  polyatomic  heats  of  formation.  Because  the  parameters  are  fitted  to  heats 
of  formation  of  molecular  species,  the  B  AC  contains  corrections  to  the  zero  point  vibrational 
energy  of  the  molecules  as  well  as  to  the  absolute  energies.  The  values  of  the  parameters  of  the 
BAC  are  therefore  dependent  upon  the  fact  that  the  geometry  and  Hessian  matrix  are  computed  at 
the  1  IF  level  as  well  as  the  fact  that  the  absolute  energies  are  computed  by  the  MP4  method.  In 
addition,  the  values  of  the  parameters  also  depend  on  the  nature  of  the  basis  sets  used  in  the  ah 
initio  calculations,  but  we  do  not  study  basis  set  behavior  in  the  present  work.  In  the  present 
paper  we  denote  this  method  HF/BAC-MP4  to  distinguish  it  from  other  possibilities  described 
below.  The  parameters  for  the  HF/BAC-MP4  method  used  in  this  study  are  presented  in  Table  1. 

Using  the  IIF/BAC-MP4  method,  the  saddle  point  for  the  reaction  is  located  at  the  HF  level 
but  the  energy  at  the  saddle  point  is  determined  by  a  MP4  calculation  followed  by  a  BAC 
correction.  Therefore,  although  the  gradient  and  Hessian  matrix  correspond  to  the  saddle  point 
(i.e.,  the  gradient  vanishes  and  the  Hessian  has  one  negative  eigenvalue)  the  energy  along  the 
minimum  energy  path  is  not  necessarily  at  its  maximum  at  the  saddle  point.  This  inconsistency  of 
the  potential  information  makes  the  use  of  conventional  transition  state  theory'  questionable  since  it 
is  based  on  a  knowledge  of  the  potential  energy  surface  only  near  the  saddle  point,  and  the 
potential  energy  can  increase  appreciably  upon  moving  off  of  the  saddle  point  along  the  minimum 
energy  path  .  The  variational  transition  state  theory  utilizes  information  all  along  the  minimum 
energy  path  and  the  saddle  point  is  not  singled  out  as  a  special  location.  The  variations  of  both  the 
potential  and  the  frequencies  of  the  bound  modes  normal  to  the  reaction  path  are  included  in  the 
VTST  calculation  and  the  fact  that  the  gradient  vector  vanishes  at  one  location  along  the  minimum 
path  is  of  no  practical  consequence  in  the  calculations. 

The  HF/BAC-MP4  method  provides  a  practical  method  for  obtaining  information  about  the 
potential  and  variational  transition  state  theory  provides  for  a  consistent  use  of  this  information  in 
calculating  a  rate  constant.  After  the  saddle  point  is  located  at  the  HF  level,  the  minimum  energy 
path  is  located  by  following  the  path  of  steepest  decents  into  the  reactant  and  product  regions.  The 
erudieni  vector  and  Hessian  matrix  are  computed  at  the  HF  level  using  the  6-31G*  basis  and  these 
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are  used  in  determining  the  minimum  energy  path.  At  each  geometry  along  the  minimum  energy 
path  the  energy  is  obtained  by  a  MP4-SDTQ  calculation  using  the  6-3 1G**  basis  followed  bv  the 
BAC  correction.  To  insure  that  the  BAC  correction  is  a  continuous  function  of  the  reaction 
coordinate,  the  BAC  correction  is  always  computed  for  bonds  which  are  being  broken  and  made  in 
the  course  of  the  reaction  and  the  criteria  given  in  eq.  (8)  is  not  used  for  these  special  'bonds'.. 

A  major  concern  in  using  the  HF/BAC-MP4  method  is  the  accuracy  of  the  computed 
potential  energy  surface  information.  One  test  of  this  is  the  comparison  of  computed  rates  with 
experimental  ones.  Tests  of  this  nature  can  sometimes  be  inconclusive  because  cancellations  of 
errors  can  lead  to  fortuitous  agreement.  Another  test  is  to  improve  the  quality  in  the  underlying  ah 
initio  calculation  to  decrease  the  magnitude  of  the  empirical  bond  additivity  correction  and  test  the 
effect  on  the  computed  rate  constant.  As  a  first  step  in  this  direction,  we  have  reevaluated  the 
potential  information  using  the  MP2/BAC-MP4  method.  In  this  method  the  geometry  of  the 
equilibrium  structure  or  transition  state  is  determined  using  a  full  second-order  Moller-PIesset 
perturbation  theory  (MP2-SD)  calculation  including  the  effects  of  all  single  and  double  excitations 
in  the  contracted  atomic  orbital  basis  denoted  6-3 1G*.  The  matrix  of  second  derivatives  is  also 
calculated  at  the  MP2  level  at  the  same  geometry.  The  minimum  energy  path  is  located  by 
following  the  path  of  steepest  descents  using  the  gradient  vectors  and  Hessian  matrices  computed 
at  the  MP2  level.  For  each  point  along  the  minimum  energy,  the  energy  is  computed  by  the  MP4- 
SDTQ  method  using  the  6-3 1G**  basis  followed  by  a  BAC  calculation.  Because  the  geometries 
and  vibrational  frequencies  are  different  at  the  HF  and  MP2  levels,  the  BAC  correction  is  also 
different  .  The  parameters  are  refitted  to  experimental  heats  of  formation  as  in  the  HF/BAC-MP4 
method  and  the  parameters  used  in  the  MP2/BAC-MP4  calculations  are  shown  in  Table  1. 

2.3.  Anhannonicity 

In  a  variational  transition  state  theory  calculation  based  upon  a  global  potential  energy 
surface  the  effects  of  anharmonicitv  in  the  potential  energy  surface  are  included  using  the 
independent-normal-mixle  approximation  by  including  the  principal  anharmonicities  in  each  mcxle. 
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When  the  potential  energy  surface  information  is  included  as  parameters  of  the  reaction  path 
i  lamdtonian.  .inharmonic  effects  can  still  he  incorporated  by  including  third  and  fourth  derivatives 
of  the  potential  along  the  independent  normal  modes.  The  effects  of  anharmonicity  are  important 
and  must  he  included  to  obtain  accurate  estimates  of  the  rate  constants. 

The  bond  additivity  correction  is  based  upon  information  about  the  potential  energv  surface 
only  through  quadratic  terms  but  is  designed  to  fit  the  experimental  heats  of  formation  of  bound 
molecules  which  include  the  effects  of  anharmonicity.  Therefore,  the  BAC  method  has 
anharmonic  effects  built  into  it  in  an  averaged  sense,  not  on  a  mode-bv-mode  basis.  The  harmonic 
frequencies  computed  at  the  HF  level  are  known  to  be  systematically  high  and  the  BAC  method 
corrects  for  this  discrepancy  as  well  as  including  the  effects  of  anharmonicity.  Anharmonic  effects 
are  most  important  for  lowest  frequency  modes  and  for  the  reactions  studied  here  the  lowest 
frequency  mode  of  the  reaction  complex  is  a  bound  vibrational  motion  in  the  interaction  region  but 
becomes  a  hindered  rotor  in  the  asymptotic  regions.  Previous  calculations  on  this  system1  have 
found  that  at  the  IIP  saddle  point  the  hindered  rotor  partition  function  and  harmonic  oscillator 
partition  function  differ  by  only  15H.  This  is  encouraging  for  computing  the  canonical  variational 
theory  rate  constants  for  which  the  dynamical  bottleneck  is  located  near  the  HF  saddle  point  where 
the  harmonic  approximation  is  expected  to  be  valid.  However,  the  calculated  tunneling  correction 
(actor  w  ill  be  subject  to  greater  uncertainties  for  lower  temperature  where  contributions  to  the 
tunneling  factor  come  from  lower  energies  which  sample  a  more  extensive  range  of  the  reaction 
c>H>rdinate.  In  the  present  study,  all  modes  are  treated  within  the  harmonic  approximation  with 
anharmonic  effects  entering  through  the  BAC.  Limits  on  the  range  of  the  tunneling  at  different 
temperatures  will  be  reported  as  a  guide  to  the  accuracy  of  the  tunneling  corrections. 
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Several  algorithms  for  determining  the  minimum  energy  path  have  been  tested''1  and  one  of 
the  must  efficient  methods  when  numerical  accuracy  of  about  1 5'V  is  required  is  the  lowest  order 
I  aler  i  >ne-step  algorithm.  This  method  uses  only  the  information  about  the  gradient  at  the  current 
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geometry  to  predict  the  next  geometry  along  the  minimum  energy  path.  In  the  present  study  the 
minimum  energy  path  was  located  using  a  gradient  following  algorithm  suggested  by  Page  and 
Mclver41'  w  hich  utilizes  both  the  gradient  vector  and  Hessian  matrix  in  determining  the  new 
geometry  along  the  reaction  path.  The  initial  step  from  the  saddle  point  was  taken  along  the 
eigenvector  of  the  negative  eigenvalue  of  the  Hessian  matrix.  Tests  were  preformed  of  the 
convergence  of  the  ICYT/SCSAG  rate  constant  with  respect  to  the  step  size  taken  in  the  gradient 
following  algorithm  and  the  step  size  at  which  the  the  energy  and  derivative  information  was  stored 
for  use  in  interpolation  of  the  parameters  of  the  reaction  path  Hamiltonian.  Since  the  Hessian 
matrix  is  used  in  the  gradient  following  algorithm,  the  two  step  sizes  are  the  same.  Decreasing  the 
step  size  from  0. 10  to  0.05  a0  was  found  to  change  the  computed  rates  at  300  K  by  less  than  2()ff 
and  a  step  size  of  0.05  a0  is  used  in  this  work. 

The  tunneling  correction  factor  computed  using  the  SCS  AG  method  is  sensitive  to  the 
extent  of  the  reaction  path.  Extending  the  reaction  path  to  0.8  aQ  towards  the  H  +  NH3  asymptote 
and  1.2  a  ,  towards  the  H-,  +  NH-,  asymptote  was  sufficient  to  converge  the  computed  rates  to 
within  10Ct  at  300  K  and  about  25%  at  200  K.  For  the  D  +  ND3  reactions  and  its  reverse  (R4), 
the  reaction  path  was  extended  out  to  0.36  a0  towards  the  D  +  ND3  asymptote  and  0.88  a0  towards 
the  D-,  +  NDt  asymptote.  This  was  sufficient  to  converge  the  ICVT/SCSAG  rate  to  better  than 
15T  at  300  K  but  uncertainty  of  about  50%  in  the  200  K  rate  remains. 

4.  Results  and  discussion 

4. 1 .  Thermochemistry'  of  the  reactants  and  products 

As  a  guideline  to  the  type  of  accuracy  expected  from  the  HF/BAC-MP4  method  in 
computing  thermochemical  data  for  bound  molecular  species  a  comparison  is  presented  of 
experimental4 1  and  computed  heats  of  formation  and  free  energy  of  formation  of  the  reactant  and 
product  species  of  the  reactions  studied  here.  Note  that  by  definition  the  heat  of  formation  and  free 
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energy  of  formation  of  II,  and  l ),  are  zero  at  all  temperatures.  The  thennoehemieal  data  lor  the 
remaining  species  are  show n  in  Tables  II. 

The  BAG  MP4  results  are  in  excellent  agreement  with  the  JANAI:  tables  for  all  species 
except  Nil,  and  NO,.  For  Nil,  and  ND,  the  results  agree  well  at  lower  temperatures  but  the 

BAC-MP4  free  energies  of  formation  systematically  underestimate  those  from  the  JANA!*  tables  at 
higher  temperatures.  The  room  temperature  heat  of  formation  of  Nil,  as  reported  in  the  JAN  AT 

tables  is  in  good  agreement  with  recent  experimental  results  of  Gibson  ct.  al4-  (45.X±().3 
kcal/mol)  and  Hack  ct.  a! 4  (45.9  kcal/mol).  Based  on  bond  additivity  relationships.  Benson'-'' 
lias  predicted  the  room  temperature  heat  of  formation  of  NH,  to  be  46  kcal/mol  and  from  Benson's 
tables,  the  heat  of  formation  of  Nil,  can  be  estimated  at  higher  temperature  to  be  52  and  5X 

kcal/mol  at  1000  and  1500  K.  respectively.  These  are  consistently  higher  that  the  JANAF  tables 
which  report  44.3  and  44.0  kcal/mol  at  these  same  temperatures.  The  thermochemistry  of  NH, 
and  ND,  above  room  temperature  remains  uncertain  and  the  BAC-MP4  method  presents  a 
consistent  procedure  for  obtaining  reliable  estimates  at  these  temperatures. 

4.2.  II  +  NH;  NH,  and  H,  +  Nil,  ->  H  +  NH,  reaction  rates 

A  comparison  of  the  potential  energy  curve  V^Ep(s)  and  the  ground-state  adiabatic 
potential  curves  V  for  the  HI7BAC-MP4  and  MP2/BAC-MP4  methods  is  given  in  Figure  1. 
Comparison  of  the  rates  for  reaction  ( R 1 )  computed  by  the  ICVT/SCSAG  method  using  the 
potential  information  obtained  by  the  HF/BAC-MP4  and  MP2/BAC-MP4  methods  is  given  in 
Table  3  and  Figure  2.  The  potential  along  the  minimum  energy  path  and  the  adiabatic  potential 
curves  show  very  similar  shapes  for  the  two  methods,  although  the  MP2/BAC-MP4  results  are 
shifted  slightly  to  the  left.  T  he  classical  barrier  heights  are  very  similar,  17.1  kcal/mol  for 
I  IF  BAG- MP4  compared  to  16.7  kcal/mol  for  MP2/BAC-MP4.  and  the  adiabatic  barriers  are  15.5 
and.  1 5.7  kcal/mol  above  the  reactant  zero  point  energies  for  the  1 1F/BAC-MP4  and  MP2BAC-MP4 


The  rate  constants  computed  using  this  potential  information  agree  well  for  temperatures 
above  3(H)  K.  At  3(H)  K.  and  below  tunneling  is  very  important  and  the  rates  are  sensitive  to  finer 
details  of  the  adiabatic  potential  curves.  At  300  K,  the  rates  computed  using  the  two  sets  of 
potential  information  differ  by  only  50ft  but  at  200  K  the  difference  increases  to  a  factor  of  3.1. 
These  results  arc  very  encouraging  for  obtaining  accurate  rates  for  a  modest  computing  effort  at  all 
temperatures  except  those  where  very  deep  tunneling  is  important. 

In  Table  4  a  comparison  is  made  between  the  different  dynamical  methods  of  computing  the 
rate  constants  for  reaction  (Rl)  using  the  potential  information  obtained  by  the  HF/BAC-MP4 
method.  The  results  of  conventional  transition  state  theory  (denoted  TST)  and  the  improved 
canonical  variational  theory  (ICVT)  neglect  quantum  mechanical  tunneling  effects  whereas  the 
improved  canonical  variational  theory  with  small-curvature  semiclassical  adiabatic  ground-state 
transmission  coefficient  (ICVT/SCSAG)  includes  tunneling  within  the  vibrationally  and 
rotationallv  adiabatic  approximation.  As  described  earlier,  the  HF/BAC-MP4  energy  along  the 
minimum  energy  path  is  not  a  maximum  at  the  HF  saddle  point;  therefore,  it  is  not  surprising  that 
the  ICVT  results  are  much  smaller  than  the  TST  ones.  At  1000  K  the  TST  results  are  higher  by 
about  70 %  and  at  300  K  they  are  higher  by  about  a  factor  of  5.5.  Tunneling  is  very  important  for 
this  reaction  contributing  99.867c,  91.37c,  47.2%  and  25.6%  to  the  thermal  reaction  rate  at  200, 
300,  600,  and  KHK)  K,  respectively,  by  the  SCSAG  method.  A  large  enhancement  in  the  rate  is 
observed  when  tunneling  is  included  -  the  ICVT/SCSAG  rates  are  larger  than  the  ICVT  ones  by 
factors  of  44.5,  7.7,  2.4,  and  1.4  at  300,  400,  600,  and  1000  K,  respectively. 

The  reaction  rates  computed  by  the  the  ICVT/SCSAG  method  with  the  HF/MP4-BAC 
potential  information  for  reaction  (Rl)  are  compared  with  the  experimental  reaction  rates1'3  in 
Figure  3.  The  overall  agreement  is  very  good  except  at  the  lowest  temperatures  which  is  discussed 
in  more  detail  below.  Arrhenius  fits  to  the  experimental  results  give  activation  energies  of 
16.0x0.2  keal/mol  for  the  temperature  range  from  908  to  1777  K,1  17.2±0.8  keal/mol  for  the 
temperature  range  from  5(H)  to  1 140  K,2  and  1 4.6±  1 .0  keal/mol  for  the  temperature  range  from 
673  to  1003  K.3  Activation  energies  derived  from  the  computed  rates  show  an  interesting  behavior 


of  initially  increasing  with  increasing  temperature  but  turning  around  and  decreasing  at  even  higher 
temperature;  e.g..  the  computed  activation  energies  are  12.9,  14.7,  17.1,  and  14.7  keal/mol  over 
the  temperature  ranges  400-600  K,  600-1000  K,  1000-1500  K,  and  1500-2400  K  respectively.  In 
general  the  agreement  between  the  experimental  and  computed  activation  energies  is  excellent; 
however,  he  theoretical  results  systematically  underestimate  the  experimental  reaction  rate 
indicating  a  systematic  overestimate  in  the  free  energy  of  activation.  The  underestimation  of  the 
rate  may  be  attributed  to  a  barrier  height  which  is  too  large;  however,  low-ering  just  the  barrier  will 
tend  to  decrease  the  activation  energy  which  is  inconsistent  with  the  experiments.  Another 
possibility  is  that  the  energetics  are  approximately  correct  (i.e.,  the  ground-state  adiabatic  barrier  is 
about  right)  and  that  the  entropic  factors  controlled  by  the  higher  lying  vibrational  levels  are  not 
accurately  described. 

The  theoretical  rate  constants  fall  within  the  experimental  error  bars  of  Marshall  and 
bontijn-  for  all  temperatures  above  660  K.  Below-  this  temperature  the  experimental  results  show 
much  more  curvature  in  the  Arrhenius  plot  than  the  theoretical  rates  and  at  500  K  the  experimental 
rates  are  higher  than  the  theoretical  ones  by  about  a  factor  of  4.  Such  curvature  in  the  Arrhenius 
plot  is  often  attributed  to  quantum  mechanical  tunneling,  but  in  the  present  calculations  the  large 
amount  of  tunneling  necessary  to  reproduce  the  experimental  curvature  is  inconsistent  with  the 
excellent  agreement  between  theory  and  experiment  for  temperatures  above  660  K.  For  example, 
at  660  K.  the  SCSAG  tunneling  correction  factor  is  about  2  and  gives  good  agreement  with 
experiment.  At  5(X)  K  the  SCSAG  tunneling  correction  factor  is  3.6  but  would  have  to  be  about 
15  to  agree  with  experiment.  To  reproduce  this  sharp  increase  in  tunneling  over  such  a  small 
change  in  temperature  would  require  a  radically  different  shape  of  the  adiabatic  potential  which  is 
inconsistent  with  the  reliability  established  in  previous  tests  and  uses  of  the  methods. 

Table  5  presents  the  reaction  rates  computed  by  the  the  ICVT/SCSAG  method  with  the 
HF7.MP4-RAC  potential  information  for  reaction  (R2)  and  these  are  compared  with  the 
experimental  ones  of  Hack  ct.  al.-'1  and  Demissy  and  Lesclaux6  in  Figure  4.  The  agreement  with 
the  low  temperature  results  of  Demissv  and  I.esclaux  is  excellent  (differences  of  18-4()(,'<  )  but  at 


higher  temperatures  the  computed  rates  are  factors  of  2.2  to  2.7  lower  than  the  results  of  Hack  ct. 
dl.  The  discrepancy  at  the  higher  temperatures  can  be  attributed  to  uncertainties  in  the  heat  of 
formation  of  Nl  1-,.  At  lower  temperatures  the  computed  heat  of  formation  of  NH-,  agrees  well  with 

the  experimental  values  and  we  expect  less  uncertainty  in  the  computed  rate  constant. 

The  rates  of  reaction  (R2)  are  related  to  those  of  reaction  (Rl)  by  detailed  balance. 

Know  ledge  of  the  equilibrium  constant  as  a  function  of  temperature  allows  calculation  of  rates  for 
one  of  the  reactions  from  the  other.  At  lower  temperatures  the  heats  of  formation  of  all  four 
species  (two  reactants  and  two  products)  are  accurately  know'n,  thus  the  equilibrium  constant  can 
be  accurately  computed  at  these  temperatures.  Therefore,  the  good  agreement  between  the 
computed  rates  and  the  experimental  results  of  Demissy  and  Lesclaux  for  reaction  (R3)  is 
inconsistent  with  the  poorer  agreement  observed  between  theory'  and  the  experimental  results  of 
Marshall  and  Fontijn-  for  temperatures  below  66 0  K.  This  is  further  evidence  of  the  inconsistency 
of  the  low  temperature  experimental  results  for  reaction  (Rl). 

4.3.  D  +  NDj  —>  Dt  +  ND-)  and  Dt  +  ND2  — >  D  +  ND3  reaction  rates 

Table  5  also  presents  the  reaction  rates  computed  by  the  1CVT/SCSAG  method  using  the 
HF/BAC-MP4  potential  information  for  reaction  (R3)  and  (R4).  These  computed  rates  for  reaction 
(R3)  are  compared  with  experimental  one  of  Marshall  and  Fontijn^  in  Figure  5.  Once  again  the 
agreement  is  excellent,  although  the  theoretical  results  systematically  underestimate  the  reaction 
rates.  This  underestimate  is  consistent  with  our  observations  for  the  H  +  NH3  — »  H2  +  NH2 

reaction  that  some  feature  of  the  potential  energy  surface  is  giving  rise  to  a  free  energy  of  activation 
which  is  too  low.  Compared  to  reaction  (Rl ),  tunneling  is  less  important  for  the  deuterated 
reaction.  For  example,  at  300  K,  the  ICVT/SCSAG  transmission  coefficient  is  44.5  for  reaction 
( R 1 )  and  only  8.5  for  ( R2),  and  at  600  K  these  decrease  to  2.4  and  1 .6,  respectively.  As  for 
reaction  (Rl ),  the  good  agreement  between  the  computed  and  experimental  activation  energy 
indicates  that  the  reaction  energetics  are  approximately  correct  but  the  entropic  factors  controlled  by 
higher  lying  vibrational  levels  leads  to  the  underestimate  of  the  rates. 
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At  the  lower  end  of  the  experimental  temperature  range,  5S0  fx,  the  agreement  between  the 
theoretical  and  experimental  rates  for  reaction  (R3)  is  excellent.  This  is  in  contrast  to  the 
divergence  of  theory  and  experiment  for  reaction  (R 11  in  this  temperature  range.  At  6(X)  K,  the 
dominant  contribution  to  the  tunneling  occurs  at  energies  about  1.1  and  0.2  keal/mol  below'  the 
adiabatic  barrier  maximum  for  reactions  ( R 1 )  and  (R3).  respectively.  At  these  energies  the  range 
of  the  adiabatic  potential  important  in  the  tunneling  calculations  is  given  by  the  classical  turning 
points  at  s=-0.35  and  0.01  for  reaction  (R1 ),  and  s=-0.22  and  -0.07  for  reaction  (R3).  For 
reaction  (R 1 )  the  tunneling  is  deeper  (at  lower  energies)  and  extends  over  a  wider  range  of  the 
adiabatic  potential  curve.  However,  the  range  of  tunneling  for  both  reactions  at  this  temperature  is 
quite  short.  To  enhance  the  tunneling  in  reaction  (R 1 ;  to  match  experiment  would  require  an 
adiabatic  potential  which  drop  off  much  more  precipitously.  This  would  enhance  the  tunneling  in 
reaction  ( R3 )  as  well  and  lead  to  an  overestimate  of  the  rate  in  this  temperature  range.  Thus,  the 
gvxxl  agreement  near  5X0  K  seen  for  reaction  (R3)  is  inconsistent  with  the  large  underestimate  for 
reaction  (Rl )  at  the  same  temperature. 


A  Conclusions 

Rate  constants  have  been  computed  for  the  gas-phase  chemical  reactions  H  +  NH;  lb  + 
N!  1-,  and  D  +  NLT  ^  Di  +  ND-,  over  the  temperature  range  from  200  K  to  2400  K.  The  rates  are 
computed  by  variational  transition  state  theory  w  ith  semielassical  adiabatic  ground-state 
transmission  coefficients1*''--  using  limited  information  about  the  potential  energy  surface  along 
the  reaction  path.  This  type  of  information  about  the  potential  can  be  obtained  directly  from  ab 
initio  electronic  structure  calculations  of  the  energy  and  its  First  and  second  derivative  with  respect 
to  ox>rdinates.  In  the  present  calculations,  ab  initio  information  is  empirically  adjusted  by  the 
B  AC-MI’4  methixi1'1'14  to  yield  more  reliable  predictions  of  the  reaction  energetics. 

The  use  of  variational  transition  state  theorv  in  these  calculations  is  mandated  for  several 
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hides  the  factors  most  important  in  controlling  the  rate  of  chemical  reactions  and 


eurrentlv  the  most  cost  effective  method  for  obtaining  reliable  predictions  of  rates  for  a  variety  of 


( 


L 


gas-phase  reaction:  (2)  it  is  capable  of  utilizing  limited  information  about  the  potential  along  the 


minimum  energy  path  v\  ithout  requiring  a  global  potential  energy  surface;  (3)  because  of  the 


nature  of  the  semiempirical  potential  information  utilized  in  this  work,  the  dynamical  bottleneck  for 


the  reaction  does  not  occur  at  the  saddle  point  for  the  reaction  and  a  variational  procedure  to  locate 


the  optimum  dividing  surface  is  needed:  and  (4)  it  provides  a  consistent  method  for  incorporating 


quantum  mechanical  tunneling  effects  which  are  crucial  for  accurate  predictions  of  the  rates. 


especially  at  temperatures  below  600  K. 


The  BAC-MP4  method  is  a  cost  effective  means  of  obtaining  reliable  information  about  the 


potential  energy  surface  that  is  needed  as  input  into  the  variational  transition  state  theory' 


calculations.  With  current  computational  capabilities,  it  is  unpractical  to  calculate  the  necessary' 


potential  information  from  ah  initio  electronic  structure  calculations,  and  the  use  of  semiempirical 


methods  is  required.  The  accuracy  of  the  semiempirical  potential  information  has  been  critically 


tested  in  this  paper.  Improving  the  quality  of  the  underlying  ah  initio  calculation  from  a  Hartree- 


Fock  (MF)  calculation  to  a  full  second-order  Nloller-Plesset  perturbation  theory  (MP2-SD) 


;alculation  decreased  the  magnitude  of  the  empirical  bond  additivity  correction.  Comparison  of  the 


rate  constants  computed  using  the  potential  information  from  the  HF/BAC-MP4  and  MP2/BAC- 


MP4  methods  showed  little  difference  except  at  temperatures  below'  3(X)  K  where  tunneling  is  more 


sensitive  to  details  of  the  potential  information.  The  good  agreement  obtained  here  is  encouraging 


for  the  applicability  of  the  more  practical  HF/BAC-MP4  to  generate  potential  information  for  a 


variety  of  chemical  reactions. 


As  a  further  tests  of  the  methods  employed  here,  the  computed  rates  were  compared  with 


experiment  for  the  reactions  H  +  N 1 1 ^  ^  H-,  +  Nil-,  and  D  +  ND^  — »  D-,  +  ND-,.  For  the  reaction 


H  +  NIC  — >  H->  +  NIC  (R1 ),  theory'  and  experiment  are  in  excellent  agreement  for  temperature 


from  660  K  to  1770  K.  For  temperatures  below  660  K  the  experimental  rates  become  larger  than 


the  theoretical  ones  with  differences  of  about  a  factor  of  4  at  500  K.  The  computed  reaction  rates 


for  the  reaction  1C  +  NIC  — >  H  +  NH^  (R2)  are  in  excellent  agreement  with  experimental  ones  but 


the  agreement  is  not  as  gtx>d  at  higher  temperatures.  The  larger  discrepancies  at  larger  temperature 


art*  attributed  to  large  iineertainties  in  the  thermochemistry  of  the  NH->  radical.  Reaction  (R2)  is  the 
res  ei  se  of  t  R 1  i  and  the  rates  are  related  by  detailed  balance;  therefore,  the  excellent  agreement 
Ivtxxcen  theory  and  experiment  for  the  latter  reaction  at  the  lower  temperatures  is  inconsistent  with 
the  larger  discrepancies  seen  for  the  reverse  reaction.  (R1 )  in  the  same  temperature  range.  The 
computed  and  experimental  rates  for  the  deuterated  reaction  D  +  ND^  — >  D->  +  NDt  are  in  excellent 

agreement  oxer  the  entire  experimental  temperature  range  from  590  to  1220  K.  Once  again,  the 
agreement  at  the  lower  temperatures  for  this  reaction  is  inconsistent  with  the  large  discrepancies 
observed  in  the  same  temperature  range  for  reaction  (Rl). 

The  overall  good  agreement  of  the  computed  rate  constants  using  different  levels  of  theory 
for  the  potential  information  and  the  good  agreement  between  the  computed  and  experimental  rate 
constants  have  given  more  confidence  in  the  theoretical  methods  utilized  here.  The  larger 
difference  between  theory  and  experiment  for  reaction  (Rl)  at  low  temperatures  indicate  that  the 
experiments  should  he  reevaluated  in  this  range. 
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Table  1.  BAC-N1P4  parameters  for  the  H  +  NH^  reaction  system. 


HE/BAC-N1F4 

MP2/BAC-MP4 

AtH({|((kcal/mol) 

5  1 .63 

5 1 .63 

AjHj^lkcal/mol) 

1 12.33 

112.53 

EN!?4(Hartree) 

-0.498232 

-0.498232 

ENlf4(Hartree) 

-34.473256 

-54.473256 

fh 

0.000001 

0.000001 

fn 

0.2174 

0.2174 

R?{(A) 

0.36 

0.36 

R;“(A) 

0.71 

0.71 

AHH(kcal/mol) 

18.98 

18.92 

ANH(kcal/mol) 

69.66 

68.06 

aHH(^'') 

2.0 

2.0 

aNH(A-’) 

2.0 

2.0 

Rn(A) 

1.4 

1.4 

Kiss 


Table  2.  Thermochemistry  for  reactant  and  product  species  for  the  reactions  NH-  +  H 
1 1-,  and  ND;  +  D  — >  NIT  +  Dt  (Energies  are  in  kcal/mol.) 


Table  3.  Comparison  of  ICVT/SCSAG  reaction  rates  (in  units  of  cm  ■'molecule' ’s'1 )  computed 
using  two  different  sets  of  information  about  the  potential  energy  surface. _ 
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Figure  Captions 


Fie u re  1. 


Ft e ure  2. 


Fie are  3. 


f'ieure  4, 


Comparison  of  potential  energy  along  the  minimum  energy  path  [V,1FP(s)]  and 


ground-state  adiabatic  potential  |  VA(s)]  as  a  function  of  reaction  coordinate  s 


computed  using  the  HF/BAC-MP4  method  (solid  curves)  and  the  MP2/BAC-MP4 


method  (dashed  curves)  for  the  reaction  H  +  NH2  -4  11,  +  NH-,. 

Rate  constants  as  a  function  of  temperature  for  the  reaction  H  +  NH2  — »  H-,  +  NH- 


:omputed  using  the  improved  canonical  variational  theory  with  small-curvature 


semiclassical  adiabatic  ground-state  transmission  coefficient  from  potential 


information  obtained  using  the  HF/BAC-N1P4  method  (solid  curves)  and  the 


MP2/BAC-MP4  method  (dashed  curves). 


Comparison  of  computed  and  experimental  rate  constants  for  the  reaction  H  +  NFF 


H2  +  NH2.  The  experimental  results  are  from  P.  Marshall  and  A.  Fontijn,  ref. 


(2),  denoted  MF  (•);  the  Arrhenius  fit  of  J.  V.  Michael,  J.  W.  Sutherland,  and  R. 


B.  Klemm,  ref.  ( 1),  denoted  MSK  (short  dashed  curve);  and  the  Arrhenius  fit  of 


W.  Hack,  P.  Rouveirolles,  and  H.  Gg.  Wagner,  ref.  (3),  denoted  HRW  (long 


dashed  curve).  The  solid  line  is  the  result  of  the  improved  canonical  variational 


theory  with  small-curvature  semiclassical  adiabatic  ground-state  transmission 


coefficient  (ICVT/SCSAG)  from  potential  information  obtained  using  the  HF/BAC 


MP4  method. 


Comparison  of  computed  and  experimental  rate  constants  for  the  reaction  H-,  + 


Ml-,  — >  H  +  NHV  The  experimental  results  are  the  Arrhenius  fit  of  W.  Hack,  P. 


Rouveirolles.  and  H.  Gg.  Wagner,  ref.  (3),  denoted  HRW  (short  dashed  curve) 


and  the  Arrhenius  fit  of  M.  Demissv  and  R.  Lesclaux,  ref.  (6),  denoted  DL  (long 


dashed  curvet.  The  solid  line  is  the  result  of  the  improved  canonical  variational 


theorx  uith  small-curvature  semiclassical  adiabatic  ground-state  transmission 


coefficient  (ICVT/SCSAG)  from  potential  information  obtained  using  the  HF/BAC- 
V1P4  method. 

Comparison  of  computed  and  experimental  rate  constants  for  the  reaction  D  +  ND^ 
— >  D-,  +  ND-,.  The  experimental  results  are  from  P.  Marshall  and  A.  Fontijn,  ref. 
(8),  denoted  MF  (•).  The  solid  line  is  the  result  of  the  improved  canonical 
variational  theory  with  small-curvature  semiclassical  adiabatic  ground-state 
transmission  coefficient  (ICVT/SCSAG)  from  potential  information  obtained  using 
the  F1F/BAC-MP4  method. 
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Appendices  B. 

Thermochemical  kinetic  analysis  of  tunneling  and  the  incorporation 
of  tunneling  contributions  in  thermochemical  kinetics 
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and  the  Incorporation  of  Tunneling  Contributions 
in  Thermochemical  Kinetics 
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Abstract.  We  partition  the  phenomenological  enthalpy  and  entropy  of  activation  for  five  hydrogen- 
atom  transfer  reactions,  H  +  1 1-,  — >  1U  +  H,  D  +  H-,  — »  HD  +  H,  H  +  D-,  — »  HD  +  D,  O(V)  +  H-, 
— »  OH  +  H,  and  OH  +  Hi  — ►  lUO  +  H  into  quasithermodynamic  and  nonquasithermodynamic 
("substantial''  and  "nonsubstantial")  contributions.  The  latter,  which  are  not  considered  explicitly 
in  previous  thermochemical  kinetic  models,  are  very  significant.  The  present  analysis  could  serve 
as  the  start  of  a  semiempirical  data  base  for  the  inclusion  of  variational  and  tunneling  effects  in 
future  thermochemical  kinetic  models. 


1.  Introduction 


From  a  knowledge  of  the  structure  and  vibrational  frequencies  of  bound  chemical  species, 
statistical  mechanics  provides  the  prescription  for  the  direct  calculation  of  thermochemical  data  such 
as  heats  of  formation  of  stable  compounds  and  heat  release  and  equilibrium  constants  for  chemical 
reactions.1  The  structural  and  energetic  information  can  be  reduced  into  thermodynamic  quantities, 
the  enthalpy  and  entropy,  which  can  be  accurately  approximated  using  bond  additivity  and  group 
additivity  relationships.  The  quasithermodynamic  formulation  of  conventional  transition  state 
theory-’-  allows  these  ideas  to  be  extended  to  tl  e  estimation  of  reaction  rate  constants. 

In  thermochemical  kinetics  a  quasiequilibrium  is  postulated  between  reactants  and  a 

transition  state  complex  and  the  conventional  transition  state  theory  rate  constant  k*(T)  is  expressed 
in  terms  of  the  associated  quasiequilibrium  constant  K  *(T)  as 


kYl')  =  ~^(CV"  KJT) 

h  P 


In  equation  ( 1 ).  kR  is  Boltzmann's  constant.  T  is  the  temperature,  h  is  Planck's  constant.  O’  is  the 
concentration  in  the  standard  Mate,  and  An*  is  the  stoichiometric  change  in  the  number  of  moles  in 
passing  from  the  reactants  to  the  transition  state.  The  equilibrium  constant  can  be  expressed  in 
terms  simile!  to  the  usual  thermodynamic  quantities 

k'<T>  -  exp (  -  A  / i  j  /  RT)  ,2) 

.  i ) 

■A’.ere  \  /  1  ^  ;s  the  conventional  transition  state  theory  approximation  to  the  standard  state  free 

energy  ol  ucmution  at  temperature  T  and  R  is  the  gas  constant.  (  The  standard  state  for  all 
equations  .md  tables  m  this  paper  is  the  ideal  gas  state  at  1  atm.)  The  free  energy  of  activation  can 


where  AyIl  and  A^S^  are  the  enthalpy  and  entropy  of  activation,  respectively,  at  temperature 
T,  and  tire  related  to  the  equilibrium  constant  by  the  standard  thermodynamic  relationship 


A 


RT  " 


d  In  k* 
dT~  " 


(4) 


•VS°r=RT 
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R  In  k  (T) 
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For  an  ideal  gas  at  a  standard  state  of  i>0  =  1  atm.  the  concentration  of  the  standard  state  is  given  by 
R"  R T.  and.  using  eq.  tl  i,  AXS^.  and  A^ll^can  be  expressed  in  terms  of  the  reaction  rate 

constant  .  For  a  bimolecular  reaction  (An*  =  -1 )  we  obtain 
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where  i%  Avogadros  number,  and  u  is  1  for  k*(T)  in  molar  units  and  (1  for  k*(Ti  in  niolcculu: 
units, 

I n  many  cases  useful  practical  accuracy  can  be  obtained  by  estimating  A ,  S  .  ■  d  \  . !  i 

by  group  additivity  and  computing  the  rate  constant  from  eqs  (li  -  i.F  '  I  ho  ;  c 

T  •» 

appealing  when  conventional  transition  state  theory,—  '  or  variational  tr.iiiMt. on. 
a  temperature-independent  transition  state,  is  accurate.  In  this  case,  the  q...i-  • 
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activation  parameters  may  be  interpreted  in  terms  of  the  free  energy  (or  enthalpy  and  entropy) 
change  on  passing  to  a  single  transition  state  complex.  For  this  interpretation  the  transition  state  is 
treated  just  like  a  bound,  stable  molecule  except  that  one  degree  of  freedom  is  missing, 
corresponding  to  the  reaction  coordinate:  all  the  bound  internal  degrees  of  freedom  can  be  analyzed 
in  terms  of  the  structure  and  vibrational  frequencies.  In  many  cases,  however,  this  simple  picture 
breaks  down  because  the  variational  transition  state  depends  on  temperature  or  is  not  a  good 
dynamical  bottleneck,  because  of  nonequilibrium  effects,  or  because  of  quantal  effects  on  reaction- 
coordinate  motion.4''  In  such  cases,  if  the  rate  constant  is  still  represented  by  eqs.  (1)  and  (2), 
A  xG'  cannot  be  calculated  by  statistical  mechanics  from  the  properties  of  a  single  transition  state 

complex.  In  the  spirit  of  group  additivity  though,  an  attempt  can  be  made  to  represent  these  effects 
by  additional  "group"  contributions.  We  will  call  such  contributions  "nonsubstantial" 
contributions. 

The  rate  constant  can  be  expressed  in  terms  of  a  transmission  coefficient  k(T)  which 
includes  the  nonsubstantial  contributions  and  a  rate  constant  kGTST(T,s)  based  on  a  single 
temperature-independent  location  sof  the  generalized  transition  state: 

k(Ti  =  k(T)  kGTST(T,s)  (8) 

We  intnxluce  the  factorizations 


,  GIST 
k  (  1 .  s) 


-<cVn  exp(AsS^/  R)  exp  (  -  AsIl”  /  RT) 


MT.  =  exp(AnS,’./R)exp(-Ann“/RT) 


(10) 


where  AsS°  and  ASH°  are  substantial  activation  parameters  obtained  from  the  properties  of  a 

generalized  transition  state,  "a  substance",  by  expressions  analogous  to  eqs.  (4)  and  (5),  and 
AnSCj.  and  AnH°  are  nonsubstantial  activation  parameters  obtained  from 


o  nT  2  d  In  k 

AnHx=RT  ~ dT 

A  S°  =  RT  +  R  In  k(T) 

n  T  dT 

The  nonsubstantial  activation  parameters  AnS°  and 


temperature-independent  and  temperature-dependent  factors  in  the  transmission  coefficient, 
respectively.  Since  AnS°  and  A nH°  are  quasithermodynamic  parameters,  their  temperature 


dependences  are  given  by  the  usual  formulas 


AsS°=AsS°  +  (AsC°p/T)  dT 


A.H°=AJf°  +|I  AT°DdT 

S  1  5  I  1  ^  r 


(13) 

(14) 


where  A^C ”  is  the  temperature-dependent  substantial  heat  capacity  of  activation.  The  parameters 
AnS”  and  AnH°  are  purely  phenomenological;  however,  because  of  their  definition  in  eqs.  (11) 

and  (12),  their  temperature  dependences  can  also  be  obtained  from  a  single  heat  capacity  analog 
AnCp  given  by  expressions  analogous  to  eqs.  (13)  and  (14).  In  this  case  the  nonsubstantial  heat 

capacity  cannot  be  given  the  physical  interpretation  of  the  change  in  heat  capacity  of  the  molecular 
complex  as  it  proceeds  from  reactants  to  the  transition  state,  but  it  will  be  useful  for  providing  a 
compact  representation  of  the  temperature  dependence  of  both  AnS°  and  AnHT_ 

Estimates  of  the  nonsubstantial  activation  parameters  can  be  combined  with  estimates  of  the 
quasithermodynamic  activation  parameters  to  yield  more  accurate  rate  constants.  In  order  to 
estimate  the  nonsubstantial  contributions  in  the  general  case  it  is  useful  to  have  some  experience 
with  their  magnitude.  In  the  present  article  we  consider  five  reactions  for  which  variational 


transition  state  theory  calculations8'15  with  transmission  coefficients  to  account  for  quantal  effects 
on  the  reaction  coordinate  motion  provide  rate  constants  in  good  agreement  with 
experiment,15,16'27  and  we  analyze  the  theoretical  rate  constants  to  obtain  values  for  AnS° , 
AnH^,  AsS°  and  ASH°.  The  reactions  considered  are 


H  +  H-,  — >  H-,  +  H 


D  +  H-,  — »  HD  +  H 


H  +  Di  — »  HD  +  D 


0(5P)  +  H-,  ->  OH  +  H 


OH  +  -»  H^O  +  H 


In  these  cases  the  largest  contributions  to  the  nonsubstantial  activation  parameters  are  due  to 


tunneling. 


2.  Theorv 


We  consider  two  factorizations  of  the  total  variational  transition  state  theory  rate  constant. 


The  first  allows  a  quantitative  assessment  of  the  improvement  over  conventional  transition  state 


theory  obtained  by  including  the  effects  of  variationally  locating  the  transition  state  and  the  effects 
of  quantum  mechanical  tunneling.  This  factorization  is  based  upon  the  improved  canonical 
variational  theory10  (ICVT)  rate  constant  k,CVT(T)  with  a  semiclassical  transmission 
coefficient10,14,28'51  k[u[)(T)  that  accounts  for  quantal  effects  (tunneling  and  nonclassical 


reflection)  on  reaction  coordinate  motion 


k(T)  =K-Uin(T)k^v,(T) 


=  KuB.cn  Kvar(T)  k*(T) 


(16) 


where  Kvar(T)  is  defined  by 

Kvar(T)  =  kICVT(T)/k*(T) 


(17) 


Both  Kvar(T)  and  Klun(T)  are  parametrized  as  in  eq.  (10)  leading  to 


^aS°T  =  AlunS(;+AvarS0T  +  A^S0T 
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(18) 

(19) 


where  A^S”  and  A^H^  are  given  by  eqs.  (6)  and  (7)  and  represent  the  substantial  contributions 
and  A[unST,  A^H^,  AvarST,and  A  varH°  represent  nonsubstantial  contributions. 

The  second  factorization  attempts  to  maximize  the  substantial  contributions  to  the  enthalpy 
and  entropy  while  minimizing  the  nonsubstantial  contributions  by  basing  the  substantial  part  on  a 

single  temperature-independent  transition  state  located  at  the  maximum  of  the  adiabatic  ground-state 
potential  curve10  sAG  (the  variational  transition  state  at  OK)  instead  of  the  conventional  transition 

state  at  the  saddle  point.  We  write 

kIC VT(T)  =  Klhcmi(T)  kGTST(T,sAG)  (20) 


where  Klhcnn(T)  accounts  for  finite-temperature  deviations  due  to  the  temperature  dependence  of 
the  improved  canonical  variational  transition  state.  Analysis  of  kGT^T(T,sAG)  by  eq.  (9)  yields 
AsS”  and  ASH”  and  analysis  of  eqs.  (15)  and  (20)  yields 


*  C°  *  nH  .  „<> 

AnS  =  A  tM  S_  +  A  .  S_ 
n  T  tun  T  ihcrm  T 


(21) 
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In  all  cases  we  used  ab  initio  potential  energy  surfaces,  namely  the  surface  of  Liu, 
Siegbahn,  one  of  the  authors,  and  Horowitz-*2'34  for  reactions  (R1)-(R3),  a  surface  called  M2,12 
based  on  the  modified35  polarization  configuration  interaction  calculations  of  Walch  et.  «/.36"38  for 
reaction  (R4).  and  the  Walch-Dunning-Schatz-Elgersma  surface39-40  for  reaction  (R5). 

For  reaction  (R1)-(R4),  tcUin(T)  is  approximated  by  the  least-action  ground-state  (LAG) 

method.31  In  these  cases  all  geometries  along  the  minimum  energy  paths  are  collinear,  stretching 
anharmonicity  is  included  by  the  WKB  approximation41  for  the  ground  state  and  the  Morse  I 
approximation42  for  excited  states,  and  bending  anharmonicity  is  included  by  a  quadratic-quartic 
approximation.43-44  For  reaction  (R5),  tctun(T)  is  approximated  by  the  small-curvature 

semiclassical  adiabatic  ground-state  (SCSAG)  approximation. I4-28'31  For  this  reaction  all 
geometries  along  the  minimum  energy  path  are  coplanar,  anharmonicity  is  included  by  the 
independent-normal-mode  approximation6-14  by  the  Morse  III  approximation14  for  in-plane 
vibrations  and  by  a  quadratic-quartic  approximation  for  out-of-plane  bends. 

3.  Results  and  Discussion 

Tables  I-V  show  results  of  the  factorization  of  the  activation  parameters  into  contributions 
from  conventional  TST,  variational  effects,  and  tunneling  effects  for  the  five  reactions  over  a  range 
of  temperatures  from  200  to  15()0K.  Table  VI  summarizes  the  results  of  the  factorization  into 
substantial  and  nonsubstantial  contributions  at  four  selected  temperatures  for  all  five  reactions.  The 
finite-temperature  deviations  of  the  VTST  contributions  from  their  substantial  part  are  also  shown 
for  reference.  Table  VII  provides  a  compact  presentation  of  the  substantial  and  nonsubstantial 
contribution  to  the  activation  parameters  for  the  five  reactions  in  terms  of  the  heat  capacities  of 


The  nonsubstantial  contributions  (primarily  from  tunneling)  to  both  the  enthalpy  and 
entropy  of  activation  are  seen  to  be  quite  large.  For  the  systems  studied  here  the  classical  barriers 
are  fairly  high  (9.8  kcal/mol  for  H  +  H->  and  isotopic  analogs,  12.6  kcal/mol  for  O  +  II->,  and  6.1 
kcal/mol  for  OH  +  H->)  which  favors  locating  the  variational  dividing  surface  close  to  the  saddle 

point.  Therefore,  the  variational  contributions  to  the  enthalpy  and  entropy  of  activation  are  in 
general  small;  they  are  largest  for  the  reaction  with  the  smallest  barrier,  OH  +  H2,  where  they 
contribute  up  to  0.5  kcal/mol  to  AvarH°  and  3  cal  mol'^K'1  to  AvarS° .  The  high  barriers  also 


var~  x- 


lead  to  considerable  tunneling  at  low-  temperatures.  This  is  reflected  in  the  large  contributions  to 
the  enthalpies  and  entropies  of  activation  at  temperatures  from  about  600K  and  below.  Tunneling 
effectively  lowers  the  threshold  for  reaction  below  the  classical  threshold,  and  thereby  makes  a 
negative  contribution  to  the  enthalpy  of  activation.  Thus,  since  tunneling  decreases  with  increasing 
temperature,  A[unH°T  also  increases  (becomes  less  negative).  We  see  a  similar  trend  for  the 
entropy  of  activation;  A  S°  is  consistently  less  than  zero  but  decreases  in  magnitude  for 


increasing  temperature. 


Table  VI  shows  that  most  of  the  effect  of  variationally  optimizing  the  location  of  the 
transition  state  dividing  surface  at  each  temperature  can  be  obtained  by  placing  the  transition  state 
dividing  surface  at  a  single  temperature-independent  location  -  the  maximum  of  the  ground  state 
adiabatic  curve  (s=s ,  ).  A  measure  of  the  importance  of  placing  the  dividing  surface  at  a  location 
different  than  sA„G  is  provided  by  ^thcrm^T  anL*  ^  thermal  w^ic^  are  b01*1  small  -  these 
quantities  vary  from  -0.1  to  0.0  cal  mor’K"1  and  from  -0.1  to  0.1  kcal/mol,  respectively,  for  all 
five  reactions  from  200  to  1500K.  The  effect  of  the  variational  contributions  in  Tables  I-V  are 
approximately  combined  with  the  TST  contributions  to  give  the  substantial  contributions  in  Table 
VI  (e.g..  ASH”  =  A  *H"  4-  A  varH”  )  and  the  nonsubstantial  contributions  are  dominated  by  the 

tunneling  contributions. 

Another  interesting  trend  seen  in  Tables  I-V  is  that  A  ,H°r  monotonically  decreases  with 

increasing  temperature  for  all  reactions.  From  eq.  (6)  w'e  see  that  if  lnk*(T)  varies  linearly  with 
1/T  then  A  ^n!j.  will  decrease  with  temperature  from  the  -2RT  term.  At  the  low  temperature  the 


r 


V. 
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decrease  is  approximately  equal  to  -2RT  but  at  higher  temperatures  the  curvature  of  Ink* (T) 
becomes  larger  and  A  ,11^,  decreases  less  rapidly.  Because  the  variational  contribution  are  small 

this  trend  is  also  reflected  in  the  substantial  enthalpy  of  activation  seen  in  Table  VI.  Also,  as  noted 
above,  the  nonsubstantial  contributions  have  the  same  trends  as  the  tunneling  contributions. 
Because  tunneling  is  so  important  at  the  lower  temperatures,  the  total  enthalpy  of  activation 
(substantial  plus  nonsubstantial)  increases  with  increasing  temperature  for  all  five  reactions  at  the 
lowest  few  temperatures.  For  all  systems  except  OH  +  Hi  this  trend  reverses  at  higher  T  so  that  the 

total  enthalpy  of  activation  does  not  have  a  monotonic  dependence  on  the  temperature. 

These  trends  are  also  exhibited  in  the  heat  capacities  shown  in  Table  VII.  We  see  that  for 
all  five  reactions  the  substantial  contributions  to  the  heat  capacity  of  activation  are  negative  and  the 
nonsubstantial  contributions  are  positive.  At  low  temperatures  (T<40()  K)  the  magnitude  of  the 


nonsubstantial  contributions  to  AC  p  are  larger  than  those  of  the  substantial  contributions,  leading 


to  a  net  increase  in  both  AaS  _  and  AaH  for  increasing  temperature.  For  all  reactions  except  the 


OH  +  Hi  reaction  the  magnitude  of  the  substantial  contributions  to  AC  becomes  larger  at 


temperatures  above  600  K.  and  both  AaS  ^  and  AaH  T  decrease  with  further  increases  in  T. 
It  is  typical  to  write  the  rate  constant  in  Arrhenius  form 


k(T)=  At  exp  [  -  Ka( T)  /  k  RTj  (23) 

and  express  the  Arrhenius  parameters  in  terms  of  the  enthalpy  and  entropy  of  activation.  From  the 


folman  definition 


4547 


at  the  activation  enertiv  we  obtain 


!Vk,(T>=  AnH‘;+  AnIil’.+  2RT 


(24) 


and  using  eqs.  ( X  H  10).  (23 ).  and  ( 24 )  we  obtain 
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For  all  five  reactions  studied  here,  the  substantial  contributions  to  the  activation  energy  are  nearly 
constant  at  low  temperature,  changing  by  less  than  0.2  keal/mol  from  200  to  400  K.  I  lowever, 
because  of  the  importance  of  tunneling  at  low  temperature  the  total  activation  energies  are  greatly 
decreased  at  these  lower  temperature  -  by  as  much  as  5.7  keal/mol  at  200  K  for  the  O  +  Hi 


reaction  -  and  the  activation  energy  changes  significantly  over  the  entire  temperature  range.  For 
example,  the  activation  energy  for  the  O  +  Hi  reaction  is  4.8,  6.4,  8.2,  9.9,  11.6,  and  13.3 

kcal/mol  at  200,  300,  400,  600,  100,  and  1500  K.  Similarly,  the  change  in  the  entropy  of 
activation  with  temperature  and  the  explicit  temperature  dependence  in  eq.  (25)  lead  to  a  change  in 
A  t-  for  this  reaction  of  4  orders  of  magnitude  from  200  to  1500  K. 

Using  estimates  of  the  contributions  to  the  conventional  TST  entropy  of  activation  based  on 
group  additivity  relationships,  Benson-'  has  approximated  the  Arrhenius  A  factor  for  the  H  +  Di 
reaction;  an  estimate  for  A^S1’  of -18.7  cal  mol^K'1  yielded  A.(|()  =  10'9  xcm'molecule'1s'1. 
The  estimate  of  A*S1^  is  slightly  higher  than  our  computed  value  of  -20.6  cal  mor'K'1  but 

neglects  the  contribution  of  -3.2  cal  tnol'1K.'1  from  tunneling.  Combining  the  substantial  and 
nonsubstantial  contributions  yields  our  computed  value  of  A  r  =  10‘l()-9cm^molecule'ls‘1.  Our 

results  show  that  tunneling  contributions  can  change  the  Arrhenius  A  factor  by  over  an  order  of 
magnitude  at  low  temperatures.  Benson's  estimate  of  the  Arrhenius  A  factor  at  300K  is  in  better 
agreement  with  our  computed  value  at  1000K  ( lO'^cm-Snolecule'V1 ),  where  tunneling 
contributions  are  less  than  3()9f. 

4.  Conclusions 

A  new  partitioning  of  the  phenomenological  enthalpy  and  entropy  ot  activation  into 
quasithermodvnamic  and  nonquasithermodynamic  contributions  has  been  presented.  The  former 
are  obtained  from  the  properties  of  a  single  temperature-independent  transition  state  (a  "substance") 
and  are  called  substantial  contributions.  The  nonquasithermodynamic  contributions  are  obtained 
from  more  global  properties  of  the  potential  energy  surface  (including  the  dependence  of  the 
variational  transition  state  on  temperature  and  quantum  mechanical  tunneling  effects)  and  are  called 
nonsubstantial  contributions.  An  analysis  of  these  contributions  has  been  performed  on  five 
reactions  for  which  the  potential  energy  surface  is  known  to  give  computed  rate  constants  in 
excellent  agreement  with  experiment.  This  analysis  shows  that  the  nonsubstantial  contributions 


can  be  very  significant;  for  example,  for  the  OH  +  H->  reaction  the  substantial  and  nonsubstantial 

contributions  to  the  enthalpy  of  activation  at  300  K  are  4.8  and  -3.8  kcal/mol,  respectively,  and  for 
the  O  +  I  U  reaction  the  substantial  and  nonsubstantial  contributions  to  the  entropy  of  activation  at 

300  K  are  -20.9  and  -8.9  cal  mol'K"1,  respectively,  for  a  standard  state  of  1  atm. 

All  the  reactions  studied  here  are  hydrogen  (or  deuterium)  atom  transfer  reactions  with 
fairly  high  barriers.  For  these  reactions  the  nonsubstantial  contributions  are  predominantly  from 
quantum  mechanical  tunneling  and  the  effect  of  variationally  optimizing  the  location  of  the 
transition  state  dividing  surface  is  small.  For  reactions  with  smaller  barriers  we  expect  the  effect  of 
quantum  mechanical  tunneling  to  be  smaller,  but  for  those  reactions  the  effects  of  variationally 
locating  the  dividing  surface  will  become  more  important  especially  at  higher  temperatures.4^’49  It 
is  therefore  expected  that  the  nonsubstantial  contributions  will  be  significant  for  a  wide  range  of 
gas- phase  chemical  reactions. 

The  temperature  dependence  of  the  quasithermodynamic  (substantial)  contributions  to  both 
the  enthalpy  and  entropy  of  activation  can  be  expressed  in  terms  of  a  single  temperature-dependent 
substantial  heat  capacity  of  activation.  Because  of  the  definitions  of  the  nonquasithermodynamic 
(nonsubstantial)  contributions  they  can  also  be  expressed  in  terms  of  a  single  temperature- 
dependent  nonsubstantial  heat  capacity  of  activation.  This  allows  the  enthalpy  and  entropy  of 
activation  to  be  compactly  tabulated  in  terms  of  their  values  at  a  single  temperature  (e.g.,  room 
temperature)  and  the  temperature  dependence  given  by  the  heat  capacity  of  activation  at  several 
temperatures.  This  provides  the  basis  for  establishing  a  semiempirical  data  base  with  includes 
important  variational  anil  tunneling  effects  in  a  thermochemical  kinetic  model. 
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.  V 


Table  I.  Activation  parameters  for  the  reaction  II  +  H-,  — »  H-,  +  H  on  the  LSTH  potential  energy 
surface.  (Rnthalpies  in  units  of  kcal  mol'1,  entropies  in  units  of  cal  mol  ’K'1,  standard  state  is  1 
atm.) 
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Table  111.  Activation  parameters  for  the  reaction  H  +  D-,  — »  HD  +  D  on  the  LSTH  potential  energs 
surface.  (Enthalpies  in  units  of  kcal  mol1,  entropies  in  units  of  cal  mol'1  K'1,  standard  state  is  1 
atm.) 
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Table  V.  Activation  parameters  for  the  reaction  OH  +  Ht  — »  H-,0  +  H  on  the  LSTH  potential 
enerc;.  surface,  tlaithalpies  in  units  of  kcal  mol'1,  entropies  in  units  of  cal  mol'1  K'1,  standard 

state  is  1  atm  1 
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Tabic  VII.  Activation  parameters  for  five  reactions.  (Knthalpies  in  units  of  kcal  mol'1,  entropies  and  heat  capacities  in  units  of  cal  mol'1 
standard  state  is  1  atm.) 


